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Abstract. Reynolds' lubrication approximation is used extensively to study flows between 
moving machine parts, in narrow channels, and in thin films. The solution of Reynolds' equation 
may be thought of as the zeroth order term in an expansion of the solution of the Stokes equations in 
powers of the aspect ratio e of the domain. In this paper, we show how to compute the terms in this 
expansion to arbitrary order on a two-dimensional, x-periodic domain and derive rigorous, a priori 
error bounds for the difference between the exact solution and the truncated expansion solution. 
Unlike previous studies of this sort, the constants in our error bounds either are independent of the 
function h{x) describing the geometry or depend on h and its derivatives in an explicit, intuitive 
way. Specifically, if the expansion is truncated at order 2k, the error is 0(£^*^^), and h enters into 
the error bound only through its first and third inverse moments h{x)~"^ dx, m = 1,3, and via 
the max norms W-^h'-^dihW^, l<e<2k + 2. We validate our estimates by comparing with finite 
element solutions and present numerical evidence that suggests that even when h is real analytic and 
periodic, the expansion solution forms an asymptotic series rather than a convergent series. 
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1. Introduction. Reynolds' lubrication equation [551 [501 [IHl [H] is used exten- 
sively in engineering applications to study flows between moving machine parts, e.g., 
in journal bearings or computer disk drives. It is also used in microfluid and bio- 
fluid mechanics to model creeping flows through narrow channels and in thin films. 
Although there is a vast literature (including several textbooks) on viscous flows in 
thin geometries, the equations are normally derived either directly from physical ar- 
guments |16| or using formal asymptotic arguments [12| . This is acceptable in most 
circumstances as the original equations (Stokes or Navier-Stokes) have also been de- 
rived from physical considerations, and by now the lubrication equations have been 
used frequently enough that one can draw on experience and intuition to determine 
whether they will work well for a given problem. 

On the other hand, as soon as the geometry of interest develops (or approaches) 
a singularity, or if we wish to compute several terms in the asymptotic expansion of 
the solution in powers of the aspect ratio e, we rapidly leave the space of problems for 
which wc can use experience as a guide; thus, it would be helpful to have a rigorous 
proof of convergence to serve as a guide to identify the features of the geometry that 
could potentially invalidate the approximation. For example, in |25j . Wilkening and 
Hosoi used lubrication theory to study the optimal wave shapes that an animal such 
as a gastropod should use as it propagates ripples along its muscular foot to crawl 
over a thin layer of viscous fluid. In certain limits of this constrained optimization 
problem, the optimal wave shape develops a kink or cusp in the vicinity of the region 
closest to the substrate, and there is a competing mechanism controlling the size of 
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the modeling error (singularity formation versus nearness to the substrate) . We found 
that shape optimization within (zeroth order) lubrication theory drives the geometry 
out of the realm of applicability of the lubrication model; however, by computing 
higher order corrections and monitoring the errors (using the results of this paper), 
we learned that cusp-like singularities are appropriately penalized by the full Stokes 
equations, yielding nonsingular optimal solutions; see [25] for further details. 

1.1. Previous work. In most of the following papers, the Stokes or Navier- 
Stokes equations are solved in a domain fig bounded below by a flat substrate and 
above by a curved boundary y = eh{x) in two dimensions, or 2; = eh{x,y) in three 
dimensions, where £ is a small parameter and the function h is fixed. These solutions 
are then compared to the solution of Reynolds' equation (or to a truncated expansion 
solution of the Stokes or Navier-Stokes equations), and the error is shown to converge 
to zero in the limit as £ — >■ 0. 

In 1983, Cimatti [5] used a stream function formulation to compare the solution 
of Reynolds' equation to that of the Stokes equation in two dimensions. The key idea 
of the proof, which all subsequent studies (including this one) also use, is that the 
Poincare-Friedrichs inequality holds uniformly as £ — > for the rescaled biharmonic 
equation (where the domain D, — fle=i is held fixed and the equations contain the 
small parameter). Cimatti assumes h has four weak derivatives (whereas, we require 
only h G C^'^) and shows that for any compact set K C ft, 

(1.1) ||£u-M||i2(o) <Ce, ma.x{\\e^px - p^\\L2f^K),\\i^'^Py\\L^K)) < C'£^/^ 

where u is the x-component of velocity, p is the pressure, a bar denotes the solution of 
Reynolds' equation, and C is independent of e but depends on h in the first inequality 
and on h and K in the second. The scaling here in not standard: he imposes the 
boundary condition eu{x, 0) = u{x, 0) = const, which accounts for the extra factor 
of £ in each of the left-hand sides of (|l.ll) . There are a few problems with Cimatti's 
analysis, notably the dependence of C on L (the "arbitrary cutoff" used to make the 
unbounded domain bounded) and the fact that some of his arguments seem to require 
£ to be small in comparison to C~^; however, his basic approach is interesting and 
inspired much of the work that followed in this subject. 

In 1986, Bayada and Chambat [3^ generalized Cimatti's work to three dimensions. 
They analyze the Stokes equations directly rather than using a stream function for- 
mulation, assume less regularity of h (apparently only h £ C^), and state their results 
in terms of limits (i.e., the quantities , edxuf, dyuf, and p^ in the solution of the 
Stokes equations converge in to the corresponding quantities in the solution of 
Reynolds equations as £ 0); hence, they do not give rates of convergence. In a later 
paper they also studied the asymptotics of the solution at a junction between a 
three-dimensional Stokes flow and a thin film flow. 

In 1990, Nazarov [18] generalized previous work to the case of the Navier-Stokes 
equations and also showed how to treat higher order corrections in an asymptotic 
expansion in the small parameter e. He proved that if h{x,y) is smooth, then there 
is a constant C depending on h, N, and the boundary conditions such that 

(1.2) ||u - u^ll^, + \\p - £- v^iL. < c£^-l/^ 

where (u,p) is the solution of the Navier-Stokes equations, and p^ are the terms 
of the asymptotic expansion truncated at the iVth order (including a boundary layer 
expansion near the lateral edges of the thin domain) , and the norms are taken on the 
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thin domain fl^ (rather than the rescaled domain f2). As a corohary, if the expansion 
is computed with "superfluous" terms that are afterwards treated as remainders, he 
obtains the optimal estimate 



(1.3) 



.N\ 
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.1/2 
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Nazarov's paper is concise to the point of being impenetrable at times. We interpret 
""1/2^ ~ {9x,dy, e^/'^dz), but this symbol was not defined and may actually be a vari- 
able coefficient operator that incorporates the boundary conditions in its definition. 
We are also unsure of the definition of p andp^, as we would have expected p — 
to appear together. 

In a later paper [T^ , Nazarov studies the asymptotics of the solution of the Stokes 
equations in a domain in which two smooth surfaces meet at a point. This problem is 
also studied in a recent paper of Ciuperca, Hafidi, and Jai in [9] . This singular limit is 
interesting in that deriving even the first correction to the zeroth order approximation 
in the asymptotic expansion remains an open problem. 

Assemien, Bayada, and Chambat [5] have studied the important question of the 
effect of inertia on the asymptotic behavior of a thin film flow, which can in many 
cases be significant, requiring that the Navier-Stokes equations be used in place of 
the Stokes equations as the underlying model for the asymptotic expansion. We also 
mention that there is a large body of literature on the long-time behavior of solutions 
of the Navier-Stokes equations on thin domains; see, e.g., pT[ IT7] . 

In 2000, Duvnjak and Marusic-Paloka [IT] showed how to rigorously analyze the 
lubrication approximation of the Navier-Stokes equations for a slipper bearing in a 
circular geometry. The focus of their paper is on formulating the problem in cylindrical 
coordinates and showing how to adapt the zeroth order case of Nazarov's proof to 
handle the change of variables. Elrod's pioneering 1960 paper [T^] is also concerned 
with the (formal) relationship between the Navier-Stokes equations and Reynolds' 
equation for this geometry. 

1.2. Motivation and summary. None of the studies described above shows 
how the constant C bounding the error depends on the function h{x) describing the 
geometry. This is because most theorems of analysis give constants that depend on 
the domain $7, which is usually fixed. But in our case, the data h(x) of the problem 
actually specifies the domain; therefore, to obtain bounds that are independent of h, 
one must avoid or modify standard arguments for flattening the boundary, etc., so as 
not to lose track of h{x) in the analysis. Moreover, arguments based on the closed 
graph theorem or Rellich's compactness theorem must be avoided entirely, as these 
also depend on the geometry. This forces us to look for new ways to analyze old 
problems using tools that furnish explicit constants. 

In this paper, we consider only the two-dimensional, periodic Stokes equations 
with a specific choice of boundary conditions, but we derive error estimates that 
depend on h in an explicit, intuitive way. Our main result is summarized in Theo- 
rem l4.1ll which may be stated as follows: Let T = [0, l]p be the periodic unit interval. 
Uk>0,he C'^''+^'^{T), 0<ho< h{x) < 1 for xGT, and e < ro/3 (defined below), 
then the error in truncating the expansion of the stream function, velocity, vorticity, 
and pressure (in appropriate e-weighted Sobolev norms) at order 2k (keeping in mind 
that only even powers of e appear in these expansions) is bounded by 



(1.4) 
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where Vq and Vi are prescribed tangential velocities on the lower and upper boundaries 
of the domain 



and pk, Ok are constants independent of h. The bound on pressure has another term 
involving ho; see (|4.110|) below. 

The constants in (|1.4p have been divided into two types: those that are (1) given 
in the problem statement or easily computable from h; or (2) difficult to compute but 
universal (independent of h). We show how to compute the constants in the latter 
category (p^ and 9^) in section |4l see Table |4^ The constants in the former category 
(rfe and Im) help us understand the competing mechanism of singularity formation 
versus proximity to the substrate: the curvature and higher derivatives are allowed to 
diverge as long as the gap size simultaneously approaches zero in such a way that the 
homogeneous products ■j^h^~^d^h remain uniformly bounded. Although the factors 
y/li and ylzjli in ()1.4p also diverge in this limit, the norm of the exact solution 
diverges at a similar rate — so the relative error in the expansion solution truncated 
at order 2k is 0(e^'^+^), with pkVk serving as an effective radius of convergence. 

The framework we have chosen for this paper is intended to be general enough 
to cover many interesting applications (such as a crawling gastropod j25j or an "un- 
wrapped" slipper bearing) but simple enough to obtain explicit detailed estimates 
that reveal the dependence of the error on the geometry h{x). We also wanted to 
determine whether there might exist geometries for which the asymptotic expansion 
yields a convergent series. Although we do not have a rigorous proof, the answer 
appears to be negative even for the simplest case of a real analytic function such as 
h{x) = I + |sin27ra;, for which the in (jl.Sp are bounded away from zero. It is 
hoped that this work will serve as a useful first step toward obtaining similar error 
estimates for three-dimensional problems that include more general boundary condi- 
tions, incorporate end effects near the lateral edges of the domain (which we avoid by 
studying the periodic case), and include the effect of inertia or viscoelasticity. 

1.3. Outline. In scction[21 we derive Reynolds' lubrication approximation in its 
primitive and stream function formulations. In section [31 we show how to compute 
successive terms in an asymptotic expansion of the stream function. In section 13.21 
we prove a structure theorem describing the dependence of these terms on h{x) and 
its derivatives. 

In section |4l we formulate the problem weakly and analyze the truncation error 
equation using weighted Sobolev spaces and a uniform Poincare-Friedrichs argument. 
The first challenge is to find the right weighted norms on the lower and upper bound- 
aries (equivalent to H^^'^iVo) and H^^^{Ti) for fixed e) to yield manageable error 
estimates in terms of h when we change variables to straighten out the boundaries. 
In section 14.41 we reduce the problem of bounding the truncation errors to that of 
bounding the second and fourth derivatives of the two highest order terms retained 
in the asymptotic expansion, namely, HV'i'x^'' ||o and ||/i^V'2:s^s^''llo- We then use the 
structure theorem of section [3.21 to compute these norms in order to obtain the con- 
stants Pk and Ok in (|1.4p for < A; < 25. In section 14.51 '^e show how to compute 
the error in velocity, vorticity, and pressure from that of the stream function. This 
requires that we determine how the Babuska-Brezzi inf-sup constant /3 depends on 
h{x)\ see [H]. 
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ERROR ESTIMATES FOR REYNOLDS' APPROXIMATION 



5 



In section [SJ we validate our results by comparing to "exact" solutions (computed 
using finite elements) for a geometry typical of engineering applications. The result of 
this comparison is that the effective radius of convergence rkPk is within a factor of 3 of 
optimal for fc = 5, fc = 10, and perhaps all fc > 5. These calculations also suggest that 
even when h{x) is real analytic, the expansion solution is an asymptotic series rather 
than a convergent series. This is because the constants pk converge to zero as fc — > oo. 
Fortunately, pk initially increases and does not become smaller than po = 0.197 until 
2k — 26, which is already outside of the practical range of fc. Finally, in Appendix El 
we present our numerical algorithm for computing the expansion solutions, which can 
be performed symbolically using a computer algebra system such as Mathematica or 
in floating point arithmetic, e.g., in C^~^ . 

2. Reynolds' approximation. Consider the Stokes equations on a periodic 
domain of width W bounded below by a flat wall moving with constant speed Vb and 
above by an inextensible sheet moving with constant speed Vi along a fixed curve 
ri,6 = {{x, h{x)) : < X < W}; see Figure I^TTl A bar is used to distinguish a physical 
variable from its dimensionless counterpart. We nondimensionalize the variables by 
choosing a characteristic speed U and height H for the problem, and set x = Wx, 
y = Hy, h{x) = Hh{x), Vi = UVi, {u,v) = u = {Uu,U^v), and p = p^p. 
The stream function ip, flux Q, and vorticity uj introduced below satisfy ip = UH^p, 
Q = UHQ, and w = fw. 

We have in mind a situation where the aspect ratio e — H /W of the physical 
domain is small. By scaling the x- and y-axes differently, we map the problem onto 
a nicer geometry, which introduces terms in the equations that vanish in the singular 
limit e -H- 0. Specifically, we wish to find x-periodic functions u,v,p defined on the 
rescaled domain 



(2.1) n = {{x, y) : 0<a:<l, 0<y< h{x)} 
such that 

(2.2) = e^U^a: + Uyy, Py = V ^ X + Vyy , Vy = "Uj; {ill Q) 

subject to periodic boundary conditions on the left and right sides of Q and 

(2.3) (u,^) |ro = (50,0), {u,v)\j,^ = {gi,hxgi) 




-p.Au + Vp = 

V ■ u = 



h{x) 



u = Voei 



x = Q 



-62 



x = W 




- + = 
-£- Aef + = 
V ■ u = 



U = VqBi 



X = 



h{x) 

-i- 



■ ei 



-62 



X = 1 



Fig. 2.1. Geometry commonly encountered in lubrication-type problems. Left: Physical coor- 
dinate system. Right: Dimensionless coordinate system ('A^ = e^9^ -\-dy). 



6 



JON WILKENING 



on the bottom and top boundaries. Here 
(2.4) goix) = Va, gi{x) ^Vi[l + e'h'{x) 



i.e., gi{x) = Vi cos6'(x), where 6 = arctan(£ft,j;) is the angle of the curve h{x) relative 
to the horizontal. Reynolds' lubrication approximation is obtained by setting £ = 
in the equations and solving 

(2.5) Px^Uyy, = 0, Vy^-Ux, u|ro = (V'o;0), u\y^ ^ {l;K)Vi. 

If we write in the form L(u;p) = (0;0;0), where L = L^°'> + £2^(2) + £42^(4) 
given by 

f-di dA (-di o\ /o o^ 

(2.6) L=\ ay + £M -dl + £M -~dl 

\ dx dy Q ) V 0/ yo Oy 

then (|2.5[) is just the zeroth order system L("'(u;p) = (0;0;0) with zeroth order 
boundary conditions (expanding go a-nd gi in (12. 3p in powers of e). The equation for 
V decouples from the others, and we find that p is independent of y and 

Integrating from to /i and solving for we obtain 

(2.8) /^(^o + T^i)-^Q, 

where Q = jj' u{x,y) dy is the volume flux through any cross section of the fluid. 
(Q is constant since V • u = and u is tangent to Fg and Fi). Since p is periodic, 
J Px dx = 0, and we find that 

(2.9) Q = -^^^^1 (im = 1^ h{xr^ dx 

Substituting p. 91) and (|2.8p into (12.71) and using Vy — —Ux, v{x, 0) = 0, we obtain the 
solution 



The vertical component v of the velocity field is customarily omitted from zeroth 
order lubrication theory as w = eUv is 0(e) on the thin geometry $7^ of Figure [2?T] 

We may also derive (|2.10p using a stream function formulation of the problem. 
Our procedure for computing higher order corrections to the lubrication approxima- 
tion and our method for estimating the error of these expansion solutions are both 
done in the stream function formulation. Let us define 



(2.11) A,^e'dl+dl, 
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In our error estimates below, we will need to consider the inhomogeneous problem 
L{u;p) = {Fi;F2;0) with boundary conditions (|2.3p . i.e., 

-AeUe +Vp = F, 

(2-12) u|ro = (ffo;0), u|p = 

V • u = 0, ^ 

Since u is incompressible, there is a stream function -0 such that 

(2.13) U = \/ X tp ^ {tpy^-ll;^), \/ X Ue ^ e'^Vx - Uy ^ -Aeij. 

It follows from (|2.12l) that ip satisfies the rescaled biharmonic equation 

(2.14) A^IP = IPyyyy + 2e^ll)xxyy + S'^lp^xxx V X F , 

with periodic boundary conditions in the a;-direction and 

(2.15) I "^"M on To, I ut^^ \ ^i, 
where Q = J^^^^^ u{0,y) dy. Since p is periodic, /p^ p2;(x, 0) = 0, i.e., 

(2.16) / ijyyy{x,0) + Fi{x,0)dx = 0. 

Jo 

Conversely, suppose we are able to find a fiux Q and a classical solution -0 of (|2.14p 
and (P?T5)) such that ((TTB)) holds. Then we define u = V x and note that (PH)) 
implies V x (A^Ug + F) = 0, i.e., the integral 

, , f , . T-,\ , , / t = unit tangent vector along \ 

(2.17) p{x,y)^ / (A^Ue + F) -tds, ^, ■ ■ u^ n\ . t f 

^ ' ^ ^' J^^ ' V path 7 jommg (0, 0) to (x, y) ) 

is independent of the path 7. A canonical choice for 7 is 

(2.18) p{x, y) = / [e^u^x + Uyy + Fi] (^, 0) + / [e^v^^ + e^Vyy + F2] (x, 7?) drj. 

JQ JQ 

Condition (|2.16p is equivalent to requiring p(l,0) = p(0, 0), from which it follows 
that p{l,y) = p{0,y) for < y < h{0), since the integrand of the second integral in 
(|2.18p is periodic in x. By construction, the variables u, p satisfy (|2.12p . where the 
boundary condition on Fi follows from the fact that ip^ + h^ipy = there; hence, 
classical solutions of the rescaled biharmonic equation yield classical solutions of the 
rescaled Stokes equations and vice versa. Reynolds' approximation (|2.10p is recovered 
if F and e are set to zero in (|2.14p - (|2.16p when solving for ip and Q; see section [37T] 

3. Higher order corrections. In this section we show how to compute succes- 
sive terms in the formal expansion of the solution of the rescaled biharmonic equation 
(|2.14p in powers of £ = H/W. For this purpose, it is convenient to manipulate the 
equations assuming they are satisfied classically. Once we obtain formulas for the 
higher order approximations, we will show (in section |4]) that they satisfy a weak 
formulation of the problem that makes it possible to obtain error estimates. See [T5] 
for background on perturbation methods in partial differential equations. 
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3.1. A recursive algorithm. Matching like powers of e in the expansion 
(3.1) [d^y + 2e^dldl + e^dt] [^("^ + e^^P^^^ + eV^"' + •••]= 0, 

we obtain the recursion 

iZ-^o) = 

^VVVV ' 

^VVVV ^^xxyy^ 

(3-2) '^VVVV ~ ~'^^xxyy ~ V'iajsa; ; fc = 2, 3, 4, . . . . 

The boundary conditions (|2.15l) become 

(3.3) 5^(2'=) = (O, g^''^ , Q(2'=) , .gf '^)) , = 0, 1, 2, 3, . . . , 

where Bip = (V'|roj'0y|ro:'0|ri,'0t/|ri) and go{x), gi{x) were defined in (|2.4I) : 

Condition (|2.16p (with i^i = 0) becomes 



\2fc 



1 



(3.5) / ilj';.fj{x,0)dx = 0, k = 0,1,2, 







If F were nonzero in (|2.14p and depended on e in such a way that V x F had an 
expansion in even powers of e, we could incorporate these terms into (13.21) and p.5p 
as well; however, we will assume F = except in section 31 where we consider the 
general case only to derive error estimates for the F = case. Let us denote the 
right-hand side of 1^ by f'^^\x,y) for fc > 0. The terms Q^^''^ in and 

(13. 3p may be computed via 

(3.6) (^(^'^■),Q(^^)) =g(/(^'=), 5^0' 

where G is defined by Algorithm 13.11 in Figure 13.11 In this algorithm, we solve 
i^yyyy = / by integrating four times in the y-direction and then correct the bound- 
ary conditions with a cubic polynomial. The formula for Q in the algorithm may 
be derived from the one for tl) as follows. As '>po^yyy{x,0) = 0, the requirement that 
lo ^fyfl^i 0) dx — is equivalent to the condition 

(3.7) = 6/' -^Q + ^^o-^o^yh + 9oh + gih 

Jo h'^ 
Solving for Q and using j h^^ dx = I3 gives the result. 

The formulas {u,v) = {-ijjy, -■tp^), w = s'^v^ - Uy, = Uyy + e'^u^x, and Py = 
s^Vyy -f e'^Vxx allow us to compute the expansions of u, lu, and p in terms of ip: 

k>0. 



(3.8) 











Vx 1 




^yy 






,/,(2fe-2) i(2k) 

s-^xx ^yy ' 










,(2fc-2) ,(2fe) 
r^xxy ' T'^yyy ' 




= 0, 


„(2) ^ _V,(0) 




_W,(2fe-4) _ J,(2fc- 
^ XXX Txyy 






rpi'^^^d^ 



i-y 
Jo 


''^{x,!]) dri. 



k>l, 

k>l, 
k>2, 

k>0. 
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Algorithm 3.1. = G{f,go,gi): 

V'o = V*/ V= Volterra operator: Yf{x,y) — / f{x,ri)drj 







1 f'^ 2'^(i{x,h{x)) ^ -ipo^y^x, h{x)) + ffo + gijx) 



^ 2hJo h{xf h{xY 
y) = Vo(a;, y) + (.9o^(a;)) 7^ 

2 

+ (3(3 - 3i/;o(a;, /i(a;)) + Vo,iy(a;, h{x))h{x) - 2gah{x) - gi{x)h{x)) 
+ (-2Q + 2tl;a{x, h{x)) - ipo,y{x, h{x))h{x) + goh{x) + gi{x)h{x)) 
return {ifj, Q) 



Fig. 3.1. Algorithm to solve ipyyyy = f , Btp = {0, go,Q, 9l), Jq ^yvyi^^O) dx = 0. 

Equation p.2p implies that dxP'y^^ = dyp'x^'^ for fc > 0; hence, differentiating under 
the integral sign in p.Sp . we see that p"x^^ and p)f^^ actually are the partial derivatives 
of ^(2*;). Finally, our choice of Q*^*"^ ensures pi^\£.,Q) d£, = i/^y^'^J {x,0) dx = 
so that p(2fc) jg periodic. 

Using Algorithm O to evaluate (V'(°\g(°)) = G(0,Vb, Vi) yields 

(3.9) Q(") = ^^^, 

. ^3 



(3.10) = (Vb/i)| + (3Q(") - {2Vo + V,)h) ^ + (-2g(o) + [Vo + V,)h 



which agrees with Reynolds' approximation ()2.10p when vS'^\ p^'^^ are computed from 
-0^"^. To compute higher order terms in the expansion, we need to study the recursion 
p.6p more closely to determine how h will enter into the formulas for Q^^*'' and 



3.2. Algebraic structure of the stream function expansion. In this sec- 
tion, we show how the terms -f/'^^'^-' and Q^^'''-' in the stream function expansion de- 
pend on h. The key result of this section is that these higher order corrections have 
a structure similar to the zeroth order formulas p.9p and (I3.10p . but the coefficient 
on each ^ now belongs to a more complicated polynomial algebra in the symbols 
Vq, Vi, /i, the derivatives of /i, and certain weighted averages of the products of h 
and its derivatives. We also present a concise representation for the correction terms 
using matrices of rational numbers that are independent of any particular choice of 
shape function h. By splitting the analysis into one part that holds universally and 
another that depends on /i in a simple way, we are able to derive useful error estimates 
governing the expansion solution truncated at any order in section HI 

Let V = hxx, ■ ■ ■] denote the algebra of polynomials in h and its deriva- 

tives over the rationals. A typical element of V might be 3 -I- ^h'^hxxhxxx- the 
generators h, hx, etc., are treated as symbols rather than functions. Thus, if h{x) 
happens to equal 1 identically, the polynomials \ — h and are nonzero in V even 
though they are mapped to zero when V is (noninjectively) embedded in C°°{T), 
the space of C°° functions on the periodic interval T = [0, l]p. If h is not smooth, 
its derivatives can still be manipulated symbolically and various subspaces (involving 
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Algorithm 3.2. (basis generation) 








for fc = 0, . . . , fco 

$fc = {ill, (or dk 
/or j = 2, . . . , ko 

for k = j,...,ko 

return {$o, • • • ,*fco} 


= 1) 




j 

12 3 4 





* 


— dk + dk-j ) 


1 

k 2 
3 
4 


• * 

• • * 

• m • ^ 



Fig. 3.2. Algorithm to find a canonical basis for each space Hk the range < k < ko. 
Here t\ hx, ■ ■ ■ ,tk ^ ^h''~^d^h. 



terms with few enough derivatives) can still be embedded in actual function spaces 
such as C^{T). 

For any monomial a = C h^" h^l^ hl^,^ • • • G P with C 7^ 0, we define its superdegree 
to be the number of derivatives present: 

(3.11) sdeg(a) = ii + 2^2 + 3^3 H . 

If a G V, we define its superdegree to be the maximal superdegree of any of its terms, 
and set sdeg(O) = —00. Since Q is a field, sdeg(a/3) = sdeg(a) + sdeg(/3) for any 
a,/3 e T'. Wc say that a is homogeneous of superdegree k if each of its terms has 
superdegree k. 

Let H CV denote the subalgebra generated by the set {h^^^^d^h : fc > 1}, i.e., 

(3.12) n - hK,, h^K,,, ...}], 
and for fc > 0, let Hk C H denote the subspace 

(3.13) Hk ~ {0} U {a G H : a is homogeneous of superdegree fc}. 

Note that "Hfe is finite-dimensional for all fc, and "Hq = Q is the set of constant 
polynomials. We will denote the dimension of T-ik by 

(3.14) dk - dim(Hfc). 

Given an integer fco > 0, we can use Algorithm 13.21 in Figure W% to construct a 
canonical basis $fc = y~p\ , . . . , '^d^s fo'^ each Hk with fc in the range < fc < fco- 
For notational convenience, let tj stand for i^h^^^dlh. As the outer loop (on j) 
progresses, <i>fc contains a basis for the subspace of Hk that involves only the symbols 
ti, . . . ,tj. Let us denote these auxiliary sets by 

(3.15) $fej = {t\' . . . : ii + 2^2 + ■ • • + jij = fc}, 1 < J < fc- 

Then $fci = {t'l}, <i>kk = and <^kj = ^-fcj-i U t^^k^jj for 2 < j < k. In other 
words, consists of together with all products of the variables ti, . . . ,tj of 

superdegree fc that contain at least one power of tj. The first several $fc and dk are 
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given by 
(3.16) 



6 



$4 = 



4 ^^x^^^ ^ ^xx ^ hxhxxx ^ hxxx 



2 ' 4 ' 6 24 
(do, rfio) = {1,1,2,3,5,7,11,15,22,30,42}, dso = 627, dgo = 204226. 

We have found empirically that the first 75000 terms satisfy jlgfqrr) < ^fe < BFfT- 
In fact, we have recently learned of the Hardy-Ramanujan formula 



exp 1 7r-\/2fc/3 J 

(3.17) dk 7= as k^oo 

for the number of partitions of the integer k. Thus, rather than 13, the base is in fact 
gTr^iTs ^ 13.001954. 

We can now describe the structure of the stream function expansion in terms of 
the shape function h. In the following theorem, V\H.2k is the tensor product of Vi 
and 'H.2k^ where 

(3.18) Vi = {0} U {a e Q[Vo, Vi] : a is homogeneous of degree 1} 

is the space of rational linear combinations of Vo and ■ Recall from (12. 9p above that 
/„, = h{x)-^ dx. 

Theorem 3.3. The terms Q^^''\ tp^^^^ in the stream function expansion defined 
by the recursion p.6p and Alaorithm \3.1\ have the form 

(3.19) 

Q(2k) ^ h^(2k) ^ J2 Q(2^)6(2fe-2£) ^ ^(2fe) ^ I2^i2k) ^ ^ Q(2£)^(2fc-2£) ^ 

£=0 1=0 



where 

(3.20) a^''Hx,y) = ^ (^) T^TT > Z^^^'^H^, y) = E /^^H^) 



2fe+3 „ 2fc+3 

yn ^ 

n— 1 ^ ' n— 1 ^ ^ 



an. 



d 



(3.21) a^f^) e ^hVin2k, f3^^'^ e H2fe. 

J2 

Moreover, = ^ ^ dx and h^^^) = ^ J,' ^ dx. 

Remark 3.4. In addition to pinning down the way in which h appears in the 
formulas for the stream function expansion, this theorem allows us to represent ?/;'^2fc) 
and Q'^^fc) using matrices of rational numbers. Explicitly, p.20p and p.2ip hold iff 
there are matrices A^^fc)^ j^C^k) .^^j^jj entries in Vi and Q, respectively, with rows 
indexed from to (2fc + 3) and columns indexed from 1 to d2k^ and containing 
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only zeros in row 0, such that 



(3.22) 



where Y2k — (1, |, . . . , (f )'^'^^'^)^ and $2fe — ■ ■ ■ , v)i^JY are treated as column 

vectors. The purpose of the zeroth row is to make it easy to convert to orthogonal 
polynomials in y/h if desired. The final statement of the theorem asserts that the 
formulas for a'^^*^-' and fe*-^'^-' are also encoded in the matrices A'^'^^^ and B^'^''\ If we 
adopt Matlab notation and denote row i of A^^*^^ by :), then 



(3.23) 
where 
(3.24) 

Note that E. 



a 



{2k) _ 



iA(2fe)(3,:)i;. 



(2fe) 



(2k) ^ f ^(2k) 



,E. 



{2k) 



'-in Jo 



^2kix) 



dx. 



(2k) (2k) 1 

is the weighted average of ip) with weight function h^"^. 



For 



example, Ei^' = (1), Ei^' (7^ ^ dx, ^ dx)^ , etc.; see (j?!^ above. 

Example 3.5. We can now represent Q*^"^ and V'*'°^ in p.9p and p.lOp by 



(3.25) 



( '\ 












1 














-2 






3 


V V 




I V 




l 


-y 



The second order terms Q*^^^ and V'^^' involve these as well as 



(3.26) a(^) = i[yo(^, 



^(2) 
^2 ' 



6(2 



/ 





\ 




( 





\ 






















-8/15 


2/15 






-11/30 


7/15 




7/15 


2/15 




19/30 


-8/15 




2/3 


-2/3 






1/3 


-1/3 


V 


-3/5 


2/5 ) 




\ 


-3/5 


2/5 / 



^(2) 



E. 



(2) 



/ 


o\ 








9/5 


-2/5 


-6/5 


-2/5 


-3 


2 


\ 12/5 


-6/5^ 



For k > 2, ^(2fe) g^j^j^ ^(2fe) g^j.g both {2k + 4) x d2k matrices with rows and 1 
containing only zeros. These matrices are universal: the shape function h enters into 
the formulas only through Y2k, ^2k, and Em'^^ in (13.22^ and p.23p . In Appendix lAl 
we show how to compute ^t^*^) and _B(2fe) directly from the lower order matrices A^^^' 
and ^(2^) with < £ < k. 

Proof of Theorem \3.3\ We saw in Example 1 3 . 5 1 ab o ve that Q^^^ and ^p^^^ have the 
desired form. Suppose ko > 1, and the theorem holds for < fc < /cq- We must show 
that it is also true for k = ko. By (|3.6p . 



(3.27) Q^^''"' 



r i _9,/,(2feo-2) 



G 



/ (2feo- 



^2) 



0, 5^°^ 

,/,(2fco-4) „ (2ko) 



ko = 1, 
fco > 2. 
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We will use the second formula for both cases with the understanding that "^^ 
should be replaced by zero. The first step of Algorithm 13.11 is to compute 'ipj^''°\ 
Using the induction hypothesis, we obtain 

i=0 £=0 

The upper limit of the last sum can be replaced by fco — 1, since we interpret 13^~^^ 
as zero. We would like to rewrite this in the form 

^ /2ko+3 „\ fco-1 /2ko-2i+3 „\ 

(3.29) 4""^ = ^ E"i""H-)|;r +E0^"M E /^r°""n-)f;r • 

^ \ Ji=4 / 1=0 \ n=4 ) 

If we use the induction hypothesis and substitute p.20p into p.28p . the operator 
V^c^y annihilates a single power of y and antidifferentiates higher powers of y twice. 
Similarly, antidifferentiates all powers of y four times. Thus, for k — and 
4 < n < 2fc + 3, we should define 

(3.30) -2/."^ (a::^^''.-") (ogr'/.-^') 

" ^' n(n-l) n(n- l)(n-2)(n-3) 

with an identical formula for /Sn'^^ in terms of /3^'^_^^2 ^^'^ Pn^-i ■ The second term 
should be omitted when fc = 1 or n = 4, and is zero when n = 5. As part of the 
induction hypothesis, we may assume that p.30p and its /? version hold for 1 < /c < fco 
as well, so that each term in the sum over £ in p.28p also has the form described in 
(|3.29p . Note that for n > and any differentiable function ^p{x)^ 

(3.31) d^{h-''ip) = /i~("+i)(/i(9a; - nh^)ip. 

By Lemmas 13.61 and 13.71 below, hdx and multiplication by both map Hk to Hk+i 
for aU fc > 0. Thus 

(3.32) 

= h[hd^- {n - 2)h,][hd., -{n- 3)/iJ (/i-^ai^l"'^) G ^Vihn2k,„ 

with similar formulas for h"d^{a^^!:l'^'' h-~"+^) and /i"a^(-^(^2fco-4)^_„+4^^ 

elude that ali'^'' and /ji^'^' have the form claimed in (I3.2ip when k — kg and 4 < n < 

2fco + 3. Finally, we obtain Q^^fco) and -^f^feo) f^om V'^^'"'^ in (jX^ using AlgorithmEIH 
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They satisfy (IXTO)) and (jX^ if we set fc = fco and define a^^''^ = 0, /^^^''^ = 0, 

2fc+3 



(3.33) 



=4 
2/C+3 

2/C+3 



Tl=4 

2fc+3 



2k 
X J 



n=4 

a(2fe) = _i_ dx, and fo'^fc) ^ _i_ ^^i^ dx. As part of the induction 

hypothesis, we may assume p.33p also holds for 1 < fc < fco- The factors of n in 
(I3.33P are due to the terms ±.tpQ^y{x,h)h in the formula for -(/j^^'^'o) in Algorithm 13.11 
The terms 3Q(2fco)|l and -2Q(2fe«)|^ in the formula for are accounted for in 

p.l9p by extending the upper limit of the sum over £ from fco ~ 1 to and noting that 
/3(o)(a;,y) = 3|^ - 2|^. Thus, V^^'"'' and Q^^fco) have the desired form, and a^n^"\ 
pi^ko) hgjQj^g ^Yie appropriate spaces, as claimed. □ 

To complete this proof, we need two simple lemmas about the spaces T-L^, (which 
also serve as the foundation for our numerical algorithm described in Appendix E)) . 

Lemma 3.6. If k >0 and (f € Hk, then hx^p £ Hk+i- 

Proof. This follows easily from the definition of in p.l3p . □ 

Lemma 3.7. If k >0 and (p G Hk, then hdx^p G Hk+i- 

Proof. If fc = 0, then hd^^ = G Hk+i. Suppose fco > 1, and the result holds for 
k < ko. Let if € be a monomial. Then there is a fc G {1, . . . , fco} and a monomial 
13 G Hk„-k such that ip = {h^-^d^h)(3. But then 

(3.34) hdx'p = (fc - l)hx'p + {h'^d'^x^^h) (3 + {h''-^d^h) {hd^P). 

Evidently, all three terms belong to Hkg+i, the third due to the induction hypothesis. 
This result can now be applied term by term for any polynomial tp G Hk. D 

4. Error analysis. To estimate the error of the expansion of tp and Q through 
order 2fc, we show that the truncation error satisfies a weak form of the rescaled bihar- 
monic equation p. 141) with data (F, go, gi) of order e^^'^'^ . We also prove a uniform 
coercivity result for the family of bilinear forms involved in the weak formulation, 
which allows us to bound the truncation error in terms of the data. 

Throughout this section, we will treat VL and T = [0, l]p as manifolds by 
identifying the points 

^- (0,y)^ (1,2/), 0<v<h{{)), 
^ ' ' T : 0-1 

and adding a coordinate chart to each that "wraps around." In particular: a function 
in C^{Vl) or C^{T) is understood to have fc continuous periodic derivatives; dVL = 
Fq U Fi; dT = 0; the support of a function (j) G C^{fl) vanishes near Fq and Fi but 
not necessarily at a: = and x — I; and the Sobolev spaces Il'^{fl) and i?Q (fi) are the 
completions of C'^{fl) and C^(f^) in the || ■ \\k norm and thus contain only x-periodic 
functions with appropriate smoothness at x = 0, 1. 
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4.1. Weak formulation of the rescaled biharmonic equation. An interest- 
ing difference between the biharmonic equation and the Poisson equation is that the 
boundary conditions in the latter are completely specified in the problem statement, 
whereas one of them (the flux Q) in the former problem must be determined as part 
of the solution. The integral condition (12.161) . which uniquely determines Q, must 
also be reformulated weakly, since it involves more than two derivatives of ip. This 
can be done [Mj by slightly enlarging the space of test functions to include functions 
that are constant along Fi (rather than equal to there). To this end, we define 

(4.2) ^^{4>eH\n) : (0,5^0) |ro = (0,0), (</., 9^0)1^^ = (const, 0)} . 

For (p, tjj in iJ^(ri), we define the bilinear form 

Jn 

= a(°) (^, (j)) + e2a(2) (^^ 0) + ^4^(4) (^^ 0). 



(4.3) 



To obtain estimates that hold uniformly in e, it will be useful to work with the 
weighted norms and seminorms 



(4.4) 



U\\l= / i^^dA, / ^l + {ei^x?dA, |^|2^ = a,(^,V), 

Jn Jn 



l,e^^U\\l + m,e^ llV^I|2.e = v/|l^llo + I^IL + IV'IL- 



For fixed e, these norms are equivalent to the usual Sobolev norms in which e is set 
to 1 in these expressions. We use x to parametrize functions defined on Fq or Fi and 
define the weighted boundary norm 

(4.5) Ml/2,e= E [l + (27rte)']'^'kP, cu^ f g{x)e-^-^^'^dx. 

k=-oo 

We equip the dual spaces and H~^{Vl)^ = [HQ(n)'^]' with the weighted norms 

(4.6) ||^||-2,s= sup |(;,^)|, ||F||_i,, = sup |(F,(w,i;))|. 

\m2., = l ll"ll!,. + l|e''lll, = l 

Since > llV'ylli.e + lk'/'a;|!i,e, the linear functional = (F, V x V') on 

satisfies \\l\\-2,e < l|F||-i,j. 

Definition 4.1 (weak solutions). Suppose 

(4.7) heC^'\T), FeH-\n)^, goeH'/\To), gi E H^/^{Ti). 
We say that {-ip, Q) e H'^{^1) x M is a weak solution of ^lA^-^J^ if 

(4.8) a,(V',0) = (F,Vx<^) 
for all (f) £ and the boundary conditions 

(4.9) (0,5o,Q,5i) 

hold in the trace sense, where Btp :— (V'|ro?'0y|ro7'0|ri,'0t/|ri)- 
Proposition 4.2. Every classical solution is a weak solution. 
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Proof. We assume -0 G C"*(il) and (|2.14l) - (|2.16p hold classically; this requires 
additional regularity for F, g^, gi, of course. If we multiply (|2.14l) by a test function 
e C^{n) n ^' and use the identity x(V x v) = V x (^v) + (V x x) • v, we obtain 



(4.10) 

= 



-A^iP + V xF)dA 



dA 



/ (V X [0(A,u, + F)] + (V X 0) . [A,ue + F]) 
I 0(AeUe + F) • t ds + / [(V X (j)), ■ (V X A,iP) + (V X 0) • F] dA 

-(A,0)(Ae0) + (Vx,/)).F] dA, 



.)-(A,V)(Vx 



■ t 



where (V x 0)^ = ((py^—e'^ipx) and the curves Fq and Fi are both oriented from left 
to right as in Figure \2A\ The conditions 



(4.11) 



ay(/)|p = 0, 



0|r„ =0, 
4> Iti = const, 

ensure that the boundary terms are zero: the first boundary term is equal to 
(4.12) (0|j,J[p(l,Ml))-p(O,/i(O))]=O 

(with p as in (|2.17p . where it was shown to be periodic), and the second is zero since 
V X = on Fo and Fi. One more integration by parts gives J^{A^(j)){Ai.4>) dA = 
a^{ip,<j)), so (ITS)) holds. Since C^{Ti) n is dense in ^I^ and both sides of (|48|) are 
bounded linear functionals of </) e this formula holds for all </) g □ 

4.2. Uniform coercivity. The following two theorems are the key to obtaining 
error estimates for the expansion solutions of section [S] 

Theorem 4.3. The bilinear form ad-,-) is coercive on ^I' (uniformly in s) with 
respect to the weighted norm \\ ■ \\2.ey i-^-, there is a constant a > such that aW^j^W^ ^ < 
ae{ip, ip) for all e > 0, e 

Proof. Without loss of generality, we may assume the characteristic height H of 
the domain was chosen so that < h{x) < 1 for < a; < 1. We now use a standard 
Poincare-Friedrichs argument j5j. Suppose G C^(il) n 5*. Then 



(4.13) 
(4.14) 



iljy{x,ri)dri 



< 



h(x) 



\fPv{x,v)\'^ dri, 



IV'(x,2/)P- 



{y ~ v)'>Pvy{j:,v)dv 



h(x) 



< y / \i^yy{x,ri)\ dr]. 



Integrating over fi, 
(4.15) 

ll^llo<% 



\^PyfdA<^\^P\l 



/ M 2 ^ "max 



In ^ 12 

Repeating this argument on the derivatives of ip yields 

1 
2 



li^yyf dA < 



2 

2,£- 



(4.16) li^ll 



UyWl + Wei^X < 



'^y\h 



-aj(-0,V') 
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so that llV'llie — f|'Je(V^: "0)- Since C (fi) n 5* is dense in ^E*, we conclude that 
(12/19)||-0||i j < a^i'ip, V') for ah -0 € * as claimed. □ 

Theorem 4.4. A weak solution tp of the boundary value problem (I2.14l) - (j2.16p 
exists and is unique. Moreover, the following estimate holds: 

(4.17) 

ll^l|2.s 



19 

< Y^I|F||-i,s 



72 + 860 



fe^hJ 



-Ae" 



1/2 



50 



1/2, e 



91 



1/2, e 



In particular, if e < ^ with r^^ = max[\\h 



x\\ooi II 2 "'"'XX II oo 



), then 



(4.18) 



19 

|^||2,e< Y^I|F||-l,e + 15 



50 



1/2,6 



51 



1/2, £ 



Proof. We begin by constructing a function ■00 G -ff^(il) that satisfies the bound- 
ary conditions (|2.15p with Q = 0. First, we map the domain to the x-periodic unit 
square R = T x (0, 1) via the transformation 



(4.19) 



M^, y) = h{x)-^'^^Jo{x, h{x)y), 0<a;<l, 0<y<l. 



We include h '^/^ here to avoid powers of h^ ^ in (I4.27p . where /iq = mino<a;<i h{x). 
We require V'oCa;, 0) = 0, -00,^(2^,0) = h{x)^^^'^ ga{x), -00(2;, 1) — 0, and ^"0,^(2;, 1) — 
h{x)'^/^gi{x). To construct such a function, we define ( e C"'(M) via 



ay)- 

and set 
(4.20) 

where 
(4.21 




0, y<-l 

y + 2y'^ + y^, -I < y < 
y~2y^ + y\ Q<y<l 
0, l<y 



Vo{x,y) = I Cfc — 7^-^ hdfe 



(k) 



(k) 



) {k)^[l + {2TTkef]^'\ [cu^dk]^ l\g^,gi]{x)h{xr^'\-^^'''-dx. 

Jo 



The value and slope of C at y = and |y| > 1 ensure that tpo satisfies the desired 
boundary conditions. Assume for the moment that each dk is zero (i.e., gi = 0). Let 
5 be the strip T x M. We may use (|4.20p to define i/jq on all of S and take its Fourier 
transform 



(4.22) 



('Ao) (k, 77) 



V'o(a;,2/)e-2"('=^+''?^) dydx = 



J -00 
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Since C is antisymmetric and supported on [—1, 1], we have 
2 

(4.23) 





2 II ~ l|2 




2 




2 






= k-o < 

2,e,R II ll2,e,S 




+ 2 


■00 


+ 

l,e,S 


■00 



2 

2,e,S 



E 



[1 + (27rfce)2 + (27rr/)2 



(■0o) (fc,??) 



/OO 
[l + (27rt)2]^ C(i) 
CO 



dt 



105 



(i?7 



90 



1/2, e 



A similar argument works if we assume go = 0, but gi ^ 0. Thus, on i?, we have 



(4.24) 



<./!!? 

2,£ -V 105 



50 



1/2, e 



51 



1/2, £ 



Next we use the formula 'ilj(){x,y) = h{xf'^'^ ^'oi^: obtain 

(4.25) "00,1;;/ ^h^^'^i'0,xy - yh^^^'^Ki^o,yy + ]^h^^^'^K'4^0,v, 

^ ~ 3 ~ 



Using Lemma [4.51 below and < h{x) < 1, wc find that 



(4.26) / %dA 



-00 {x, h{x)y)^ h{x) dxdy — / ft.^V-'o 1^ i dA, 



"'0 



r.ydA< ro,y dA, 



e'ro,xdA<2 / eX.d^ + 4e^||/i^| 



2s^^Pl,ydA < 2 / 2eVo,.,rf^ + 4£'l|/^'ll 



2 / ^o,,,rf^+o / %,ydA 



e'ro.xxdA< - / eX..'^^ + 30e1|/i^||oo / s'^j^^dA 

^ JR. /• _ "'■'^ 

+ 30£2||/i2||^ / 2e''ijl,ydA + iQe^Wh'hl^W^ I 4,ldA 



+ 30£^||/i^||oo / ro^yydA + iOe^WhlWoo / i^iydA 



+ iOe^Whi^W^ / i^lydA + 30£^||/i^||oo / %dA. 

Jr. Jr 

Note that the inverse powers of h in (I4.25P are canceled when we change variables 
from 51 to i? by the factors of h that arise from the substitution y ^ hy and from the 
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Jacobian of the transformation. Collecting terms and majorizing, we obtain 



(4.27) ||^o||2,e < (5/2 + 30e^\\h,\\l + SOe^h^Wt, + SOe^hh^^Wl) 



1/2 



■00 



2,£ 



Finally, we correct V'o by a function in 4* to obtain the weak solution -0, which must 
satisfy 

(4.28) V-V-o e a,(V'-V'o,<?!>) = (^,0) (F,Vx0) -0,(^^0,0) for all € 

Since / is a bounded linear functional on ^, the Lax-Milgram theorem implies exis- 
tence and uniqueness of the solution ip of (|4.28l) and gives the error bound 

iq 

(4.29) 11^ - ^-0112,^ < a-'\\lh^, < - (||F||_i,e + Uoh,e) ■ 

Combining this with (I4.24p and (|4.27l) and using the triangle inequality gives (I4.17p . 
where we note that |(§)^(Yf ) < 72 and 30(fi)2(f§) < 860. □ 

The following lemma was used to balance the coefficients in the terms of l|4.25p 
as much as possible. 

Lemma 4.5. For any ai, . . . ,as 

(4.30) 

(ai +02+ 03)^ < 2aj + 4al + 4a^, 



/ I I \2 / 2 I -'^0 2 I -1 r 2 I 40 2 , on / 2 I 2 I 2\ I 2 

(oi H + as) < -a^ + —02 + ISog + —04 + 30 (og + ag + a^j H —a^- 



Proof. In general, given positive real numbers 71, . . . , 7„ such that 7^ ^ < 1, 

then for all a G M", we have (X^i "^j)^ ^ Si Tj'^j- This is a consequence of the 
Cauchy-Schwarz inequality: 




(4-31) (^E«. J - [T. (77^^^) I M E77^ I I E7.«: 

One readily checks that ^ + j + j = ^ and (| 

4.3. Truncation error equation. In section |3l we showed how to construct 
successive terms in the stream function expansion by solving the recursion p.2l) - p.5p . 
Theorem 13.31 guarantees that derivatives of h higher than 2k do not appear in the 
formulas for . . . , ip^'^''^; hence, if G C^^iT), these functions satisfy (l3:2t- ((33t 
in the classical sense (with k replaced by (. and running from to fc instead of to cxi). 
Thus, iih& C2'=+4, vipproa; = V'^"^ + E^V'^^^ + ' ' ' + e^^^^'^^ Satisfies 

(A 321 A^V/^fe) - c2'=+2 ( 2^,(2'=) + 7/;(2'=-2)\ , 2fe+4^(2fe) 

The truncation error ^/'ir? = V'ea;act— V'opproa; then satisfies A'^tp^rr = —^ei^approx with 
Q^g2/c+2-j boundary data. Since the right-hand side of (|4.32l) and the boundary data 
are known in terms of h, we are able to estimate the size of "fAerr' using Theorem 14.41 
above. However, to use this theorem, we need to formulate (|4.32l) weakly. 
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We begin by showing that the -0'^^^^ satisfy a weak version of the recursion 
Suppose fc > and h G C^'^(T). Let e 5* and denote the constant value of (p on Fi 
by q. We muhiply both sides of p.2p by cf) and integrate by parts using 



(4.33) 



Jn Jo 

</> ^il'L dA= [ (^,, -q f {x, Kx))K dx 



dx. 



to obtain the recursion 



,(0) 



,0) =0, 

(4.34) aW(v/2),0) =-a(2)(^(0)^^)^ 
By (13.81) . the boundary terms in (|4.33l) combine to form 

'p'i'\xM^))+pT\^M^))h. 



(4.35) 



dx — 0, 



£ — 2,3, . . . ^ k. 



0<£<k, 



when substituted into p. 21) . Other boundary terms do not arise in (|4.33p , since <f> = 
on Fq and 4>x = 4>y — ^ on Fq and Fi. 

Now suppose h e C^'^+^'^(T), i.e., has 2fc + l continuous periodic derivatives 
and dl'^^^h is Lipschitz continuous so that d^'^^^h e L°°{T). Let {'>p exact, Q exact) be 
the weak solution of 

(4.36) A^V-^O, B^ = (0,5o,Q,5i), 

with go J ffi given in (|2.4p . and define the truncation errors and approximate solutions 



(4.37) 



(2fc) _ 



ea;act ^'approx^ 



{2k) 



(2fc) ^ , (0) , ^2 , (2) 



2fc,/,(2fc) 



)(2fc)-n ^-OC^'^) 0(2'=) - + f2/-)(2) J , 2kQ(2k) 

'err exact ^ approx, ^ approx ^ ^c^/ ^ -pc^/ 



Since ip^°\ . . . , ■0'^^''' satisfy (|TM)) and ((H^)) while aeiipexact, 4>) = for every G 
we may expand ae{'>p^3-r\ <t>) in powers of e to obtain the truncation error equation 



(4.38) 

where 
(4.39) 

and 



^ _^2fc+2| 



Ffc,V X 



50.(2,? = (0,0, g(f;),£2'=+^fc 



'a(2)(0(o),0)+eV4)(0(o),0), fc = 0, 

a(2) (0(2fe),0) +£2^(4) (0,(2fc),^) +a(4) (^(2fc-2)^^^ ^ > 
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There are many functionals Fj. e H ^(51)^ that have this action on the subspace 
(4.40) V" = {Vx0 : e *} = e : u^r + Vy ^ . 

(AeUjjt,»TO2;)el Satisfies A^V'aOTroa; = e^'^+^VxFfe 



For example, F^ = e ^'^ ^[VpLpploa; 



*-approx)e\ sdijianeh aa^. ijJapprox — t V A j: 

classicahy and, using p.Sp . may be shown to have following action on HQ{n)^: 



(4.41) 



(Ffc,(u;z;)) 



(Ffc,(u;«)) 




0)) (u,-eV)-(e-Vi",^)(e^;,)dA 



fc = 0, 



dA, fc > 1. 



This choice is suboptimal because the terms e^^ip^xy^ and e^^V'iaf diverge as e — )■ 0. 
We grouped e with Vy and with Vx due to the definition (|4.6p of |jFfc||_i^E. Instead, 
we will use the following functional, which agrees with (|4.4ip on V: 



(Ffc, (u- v)) = j^^ (V4°j) {2uy - e^vx) dA, k = 0., 

(4.42) (Ffc, (w; v)) = J^^ [^jg^^) {2uy - s^v,) + /^VEt'^^M k > 1, 

(r;-y)zi,(x,,7) 



4>[u\{x,y) 







y 



h{x) 



drj. 



Note that if e and q — <j>\ri, then = (0 — q)h ^. The purpose of the 

here is to be able to include an with V'iraK^'' in the error estimates below (to 
properly consolidate terms). Another alternative to (|4.42p that would lead to similar 
estimates below is (Ffc, (u; v)) = J^-^i—s^ipx^x^Vx — '^^yy^'^^Uy] dA. 

4.4. Error estimates. Let us assume from now on that e < ^ with 7'o = 
max(||/i3;||oo, W'khhxxW^L'^)^^ ■ Then by (j4.38p and Theorem 14.41 we have 



(4.43) 



(2k) 



< e 



2k+2 



2,e 



-I 

12' 



15 



Ik 



1/2, £ 



It remains to bound the norms of F^ and 7^ in terms of h. From (j4.42p . we have 



(4.44) 



(l)[u\{x,y) 



< 



V 



h{xy 



h{x) 



\uy(x,'q)\'^ dr] , 



where the first integral is '■"^^4'' and the second is independent of y. Hence 



(4.45) 



1 ph{x) 

Jo 



(f>[u]{x,y) 



dy dx < 



12' 



From (|442|) . we then have 



(4.46) 



< 5a^ + 







VT2) 


\\Uy\\ 






VT2 


^12 



2\ 1/2 



{Mle + \\SV\\1,) 



1/2 
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where a = \\il)'i'x\\o and b = ||/i^V'ixxx^^||o- Using ^ < fa^ + -^b"^ , we find that 



(4.47) llFfelLi^ 
FinaUy, by (14.43^ . we obtain 



, 25 ^ 7 5 3, 19 „^ „ 

<W— b^<-a+-b, — Ffc -, < 4a + 6. 

4 20 - 2 5' 12 " "-^'^ - 



(4.48) 



{2k) 



< e 



2k+2 



2,e 



(2fe) 



/j2^(2fe-2) 



15 



7fc 



1/2, e 



where the fourth derivative term should be omitted when fc = 0. In section I4.5[ the 
following bound will also prove useful: 

(4.49) 



'2k) 



-2k+2 



2,e 



2k) 



< e 



2k+2 



(2k) 



/,2^(2fc-2) 

T^XXXX 



15 



Ik 



/ll/2 



1/2, e 



The truncation error in the flux expansion satisfies 



(4.50) 



hi.) g2^i2k) 



= ^^^^';Hx,h{x)) = I {h{x) - r^)^^{x,f^)dv 



for any x G T. Using estimates similar to (I4.44p and (|4.45p . we find that 



(4.51) 



Q 



{2k) 



< 



V3 



{2k) 



2,e 



Thus, we have reduced the problem of bounding the truncation errors to that of 
checking the norms of three quantities that can be computed explicitly in closed 
form. 

We begin by attacking the boundary term \\h^^^^^k\\i/2,e- This norm was defined 
in (|4.5p above. Recall that 



-2fc-2 



(4.52) 7. = ^1 ( [1 + e'hl] - { T) 

where ("J,/^) = 1 and for ^ > 1, 

(AK-x\ ^"l/2^ _ (~l/2)(~3/2)...(-[^-l/2]) _^ 1 3 5 2^-1 

^^•^''^ ^ )~ (l)(2)...(f) 2'l"Q'"^2r- 

Taking the logarithm of this product and its inverse, one may show that 



(4.54) 



1 



^/4FTT 



< 



-1/2 



< 



1 



^/MTT' 



0,1,2,. 



Since h e C'^^^^'^{T) C C^'^(r), we know 7^ is at least Lipschitz continuous on T 
and so belongs to H^{T). As a result. 



Ik 



1^1/2 



< 



1/2, e 



Ik 



(4.55) 



/1I/2 
1 



^-3/2 

= / h-^ll + £^ ( T^Kjk + h-^/^-ik,x ) dx 



< I h-Wk + -/h-^hl^l + \e^h~Wk,x dx. 
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Since £||/ix||oo < 1/3 < 1, the binomial expansion of converges uniformly 

on r = [0, l]p. As the terms in this expansion alternate in sign, the error in truncating 
the series is smaller (pointwise) than the first omitted term. Therefore, for each a; G T, 



(4.56) 



-1/2 
k + 1 



2fc+2 



< 



V3fc + 4 



2fc+2 



Since (-2^)(~y^) = {'^/J^), by differentiating (H3^ we obtain 



(4.57) jk,x 




-3/2 



fe-1 

E 

i=0 



-2fc-2 



The terms in the expansion of [f + e^li^] "^^^ also alternate in sign, so 



(4.58) |7fc.x(x)| < Vi 



-3/2 
k 



2{k + l) 
V3fc + 4 



Combining (|435)) . (|4?56| . and (|4?58| and using ^3^,^^^ 
(4.59) 



< 



fc, we find that 



Ik 



2 

l/2,e 3fc -I- 



+ - 



5 £2/1? 







4 /l3 



f 6(fc -I- i)^/i;^' 



2l,4fe ( ''■''■a^ 



'WhxWi 



2 



5 p(2fc+2) 



20, 

y 



V 2 

(2fe-|-2) 



da; 



/1//3 



/85 









max E 

(mj)e{(l,l),(3,l),(3,2)} "'^ 



(2/C+2) 



where 



(4.60) 



f 



■ dx. 



Recall that $2/c 



h{xy 

r (2k) (2k) ^ 



(2k) 



(x) 



h{xy 



■ dx. 



{hf, \hhf-^K 



E, 



(2k 



^e^-^dfh} is a 



basis for the space •H2fc defined in (|3T3)l . is the square of the 2-norm (or second 

moment) of iff^\x) with respect to the probability measure I^^h{x)~™ dx, whereas 



in p. 241) is the expected value of tpj ' {x) with respect to this measure. Using 



(I4.59P in (I4.48P gives a bound on the error caused by failing to satisfy the boundary 
conditions in the stream function expansion. It is perhaps surprising that this bound 
can be expressed in terms of three simple integrals involving h and its derivatives. 

Remark 4.6. Ii and I3 are dimensionless quantities in (|4.59p — if /i and x still 
carried dimensions of length, an extra length scale (e.g., /imax, which is currently set to 
1) would need to be included in the Sobolev norms to allow HV'Hoi IV'li.e: and IV'li.e to 
be added together; this length scale would also appear in (|4.59p to nondimensionalize 
/1//3 in the denominator 

in (|4.48l) in terms of h. It will be 



Finally, we estimate ||V'i2:'^''||o and Wk'^tl^: 



.2,,,(2fe-2)| 



useful in our analysis to define 
(4.61) Tk- 



max 

l<f<2fc+2 
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SO that forO<^<fc + l,l<j< d2i, and m = 1, 2, 3, we have 



(4.62) 



ip) '{x) 



< r] 



-21 



E, 



(11) 



< r] 



-21 



< < r-«. 

— m J — k 



If h is real analytic as well as periodic, a standard contour integral argument shows 
that there is an r > such that ||9^ft,||oo < fc! r^'' for all k > 0. Such an r serves as 
a common lower bound for in (|4.6ip for all A: > 0. If /i is a constant function, the 
results below hold with r/j = oo and r^T^ = (i.e., the lubrication approximation is 
exact) . 

Our first task will be to bound the growth of the terms Q'-^'^^ in the expansion of 
the flux. By Theorem 13.31 and Remark 13.41 there are rational matrices Aq^'^' , A^'^' , 
with rows indexed from to 2fc + 3 and columns indexed from 1 to d2k such 

that 



(4.63) 



I?, 



fc > 0, 



£=0 



where a^^*^) = i[l/o4 ^(3, :) + V^iAi'*^(3, :)]S^'^^ and fo^^fc) ^ is(2'=)(3, :)£;^'*\ See 
Example 13.51 above for a reminder of how this works. We now use (|4.62p together 
with the fact that ju • wj < ||u||i||t«||oo for v^w to conclude that 



(4.64) 

where 
(4.65) 



,(2^) 



< 



5(2^) 



< £ < fc. 



d-2k 



z = 0,l, 



,(2fe) 



=^5: i?(^'=)(3,j) 



It follows from (|4.63p that if we increase k. 



X2k) (2k) 



for fc= 1,2,3, . 



(2k) _ (2fc) 



fe-i 

E 

£=0 



(4.66) 

then 
(4.67) 



The constants n^n^^ , n^i^^ , k^o'^^ 



(21) (2k-2e) 



via the loop 



(i = 0,1), 



fc > 0. 





<-( 







for all along with the matrices A^^''^ and B^^fe)^ gee Table HTTI 

Now that the terms Q^^^) 
and |l/i^V'i^xx^''||o- Recall that 



do not depend on h and may be computed once and 

{2k), 



Now that the terms Q^^^) have been bounded, we are ready to estimate \\i^xx ||o 

2,/,(2fe-2) 



(4.68) 



k 

^(2'=)(x,y) = :^a(^'^)(x,y)+^Q( 



(2£)^(2fe-2£) 



where a'^2fe)(2; y) and j3'''^^\x,y) have the matrix representations (I3.22p . A slight 



modification of the proof of Theorem 13.31 shows that there are also matrices 



i(2fc) 
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Table 4.1 
tif'^^ before and after loop 114.6611 . 



k 


(2fe) 


before 


(2fc) 

«1 


before 


K 


(2fc) 
2 




after 


A2k) 


after 





5.00 X 


10-" 


5.00 X 


10-01 


1.00 


X 10+™ 


5.00 X 


10-" 


5.00 X 


10-" 


1 


3.00 X 


10-'" 


5.83 X 


10-" 


8.00 


X 10-" 


7.00 X 


10-" 


9.83 X 


10-" 


2 


5.30 X 


10-" 


7.05 X 


10-" 


1.73 


X 10+"" 


1.96 X 


10+00 


2.36 X 


10+00 


3 


2.72 X 


10+00 


3.73 X 


10+00 


6.74 


X 10+"" 


8.87 X 


10+00 


1.07 X 10+°^ 


4 


1.83 X 


10+01 


3.32 X 


10+" 


4.14 


X 10+"i 


5.43 X 


10+01 


7.32 X 


10+01 


5 


2.00 X 


10+02 


3.69 X 


10+02 


4.55 


X 10+"" 


5.28 X 


10+02 


7.30 X 


10+02 


6 


3.41 X 


10+03 


6.32 X 


10+03 


7.22 


X 10+"-'= 


8.00 X 


10+03 


1.13 X 


10+04 


7 


7.77 X 


10 + 04 


1.66 X 


10 + 05 


1.54 


X 10+"=^ 


1.68 X 


10 + 05 


2.63 X 


10+05 


8 


2.69 X 


10 + 06 


5.23 X 


10 + 06 


4.69 


X 10+"« 


5.31 X 


10 + 06 


7.98 X 


10+06 


9 


1.26 X 


10 + 08 


2.31 X 


10 + 08 


1.94 


X 10+""^ 


2.31 X 


10 + 08 


3.40 X 


10+08 


10 


6.51 X 


10 + 09 


1.45 X 


10+10 


9.97 


X 10+"3 


1.18 X 


10+10 


2.00 X 


10+10 



^(2fc)^ j^(2k) ^(2k)^ ^(2fe) of dimension (2fc + 4) x d2k+2 and (2fc + 4) x d2fe+4, 

respectively, such that 



fa^^^\x,y) = h{x)-' {Y2k{x,y)f Foi^^ + V^i^i 
f3if\x,y) = h{x)-^ {Y2k{x,v)f i3^^''^^2k+2{x), 



4(2*=) 



2fe+2 



h 



o^xxxxix^y) ^h{x) {Y2k{x,y)) 



■^oif 



{2k) 



2fe+4 



(4.69) P'ix%{x,y) = h{x)-^ {Y2k{x,y)f B^^^^^2k+i{x), 

where l2fe = (1; f > • • • j i'fif'^^^)'^ ■ We can achieve significantly sharper estimates of 
IIVE'^llo and \\hS 
polynomials. Let 



llV'iaf^llo and Ij/i^V'ixxx^^llo by expressing the dependence of i/' on y using orthogonal 



(4.70) 



Y2k - ( Pr 



, • ■ • , ^2fc+3 



(!) 



be the vector of shifted Legendre polynomials [T] , which satisfy 



(4.71) 

The first several are 



Pm{x)Pn{x) dx = 



2n + 1' 



m, n > 0. 



(4.72) Fo = 1, Pi = 2x - 1, P2 = 6x^ - 6a; + 1, P3 = 202;^ 



The well-known recurrence [T] 

~ 2n — 1 ~ 
(4.73) Pn{x) {2x - l)Pn-i{x) 



n-1 



-Pn-2{x), 



n > 2, 



12a: - 1. 



can be used to construct a nested family of lower triangular matrices i?2fe of dimension 
(2fc + 4) X (2fc + 4) with indices starting at zero such that 



(4.74) 

For example, 



(4.75) 



i?o — 



Y2k — P2fci^2fe, 



o\ 

12 

1-6 6 

V-1 12 -30 20/ 



_ yT n-T 



'-2k^^2k 



( 1 





1/2 


1/3 





1/2 


1/2 








1/6 


Vo 
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The entries of R2i[ are nonnegative and have unit column sums for all A; > 0. Next 
we renormalize the shifted Legendre polynomials and define 

(4.76) P„(x, y) = h{xy^/^V2^[+TP„{y/h{x)), 

Y2k - (Po, • • • , ^^2^+3)^ = h-^^^D2kY2k, D2k = diag (VT, Vs, . . . , V^k + l) 
so that (|4.69p becomes 
(4.77) 

Pif^ ^ {Y2,{x,y)YD^,^R-,^B(''^ (h{a 

yoir^ (M^)-'/'$2fcH-4(x)) , 



xr^^^'i>2k+2{x)' 



n-1 TD-T 

^2k ^2k 



Y^h'cx^''l-{Y2k{x,y) 

h'fiifl - {Y2k{x,y)Y D-,'R-,^M''^^ (Mx)-^/^$2fe+4(x)) 
Example 4.7. From ((X^ and (IXTU|) . we have 

'-4/i2 +2M 



(4.78) = 



2 

^0 



Is' 

Yl 
2 



(K) + ^1: 



h 

■ 2hhr 



ft2 



which has the form described in (|4.68p and (|4.69p with 
(4.79) 



m m A(o) 
^0 ' ^1 ' ^0 



( ' 

















\ 




















-4 


4 


-2 


2 


18 


-12 




V 6 


-4 


6 


-4 


-24 


12 





and the form described in (|4.68l) and (14. 77^ with 

/I 2 

/ 6 6 



Dq ^Rq ^ 



m j(0) n(0) 
^0 ) ^1 ) ^0 



10\/3 10\/3 

5 -2 

6^5 6V5 
I 3 -2 

V 10V7 10V7 



5 
6 

17 



-2 
6 

-8 



10\/3 IOV3 

7 -4 

3 -2 

10V7 IOV7 





-18 



-6 



10\/3 10\/3 
-18 6 

-12 6 



IOV7 IOV7 



where we recall that $2(2;) = (/i^, ^hh^x)^ ■ The matrices Aq\ A\' , B^^^ representing 
ip'xxxx are each 4x5 matrices while $4(0:) was given in p.l6p . 

To compute || V'ia;'^'' ||o and ||/i^V'ixa;x^^||o: note that each of the expressions in 
(|4.77p is of the form Y^n=Q Pn{x.,y)wn{x), where w{x) = Sz{x), S* is a constant 
matrix, z{x) — h^^^2k+2j{x), j = 1 or 2, and m = 1 or 3. For fixed x, we have 
/o'*^^^ Pm{x, y)Pn[x, y) dy = 5„ 



It follows that 
2 



(4.80) 



nh{x) / y .1 

\'S^Pn{x,y)wn{x)\ dydx= / V'w„(a 



' dx. 
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Moreover, J2n^"(^)^ — Il'"'(^)ll2 ^ ll'S'llFll^(a;)||i, where || • ||f and || • II2 are the 
Frobenius and 2-nornis of a matrix and vector, respectively. Integrating, we have 
/o \\z{x)\\ldx = I„^\\Eif+^'^\\l. Since < dj^r^*^ for < ^ < fc + 1, we define 



(4.81) 



(2k) _ 



2k+2 



K. 



(2k) _ 



2k+2 



K. 



(2k) 



\fd^k 



+4 



K. 



(2k) _ 



2k+4 



D2^R^iAr'^ 

D2k'R2lB^''' 
n-i u-T v(2fc) 

^2k ^2k 

p-1 TJ-T r>(2k) 
^2k ^2k ^ 



i = 0,1, fc > 0, 
fc > 0, 
i = 0,1, fc > 0, 
> 0, 



SO that 



(4.82) 



d2i) 



(21) 



r'xx 



(21) -21-2 



< 



h(\Vo\4''^ + \V,\Ki 



< VhK'i'r 



(2e)-2l-A 



o<e<k, 

0<i<k-l, 
< £ < fc - 1. 



From the bound (|4.67l) on IQ^^^-^I and the formula (|4.68p for ^(^'^^ in terms of a^^*^) 
and Z?'^*-'"^^), we see that after increasing k''^^\ k'^^\ Kq^''\ k'"^^^ via 



(4.83) 



we have 



^(2fe) ^ j^(2k) 



(21) j^(2k^2l) 
2 I 



1=0 
k 



(2fc) ^ ^(2fe) 



(2£) ^(2fc-2£) 
2 I 



t=0 



j = 0,1, fc > 0, 
i = 0,1, fc > 0, 



(4.84) 



'2k) 



(2k)\ -2k-2 



1^2 ,(2/c-2) 



< Vh 



(|Voi4 



(2fc-2) 



i2k-2)\ -2k-2 

I I. 1 



fc > 0, 
fc > 0. 



.(2fe) 



In each term Q^^'^-'/?^^'^ , we have used /| //1/3 < 1 (which follows from the Cauchy- 
Schwarz inequality) to majorize h/^/h by a/7i- The constants t Ki do not 

depend on h and may be computed once and for all along with the constants kI 
and the matrices Af^'^ and B^^*^); see Tables and HJl 

Finally, we combine the boundary estimate ()4.59p with the interior estimate (j47 
to bound the truncation error via ()4.48p . In terms of r^, (|4.59p gives 
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Table 4.2 
K^^^^ before and after loop 114.831 1. 



k 


"■0 


before 


^(2.) 


before 


^(2fe) 


„(2fe) 
"■0 


after 


^(2fe) 


after 





9.95 


X 


10-01 


2.17 


X 


10+00 


2.99 


X 


10+00 


2.49 


X 


10+00 


3.67 


X 


10+00 


1 


2.33 


X 


10+00 


4.70 


X 


10+00 


7.99 


X 


10+00 


8.41 


X 


lO^ 


-00 


1.16 


X 


10+01 


2 


7.42 


X 


10+00 


1.58 


X 


10+01 


2.51 


X 


10+01 


3.14 


X 


lO^ 


-01 


4.32 


X 


10+01 


3 


4.29 


X 


10+01 


8.62 


X 


10+01 


1.19 


X 


10+02 


1.62 


X 


lO^ 


-02 


2.21 


X 


10+02 


4 


4.58 


X 


10+02 


8.71 


X 


10+02 


1.03 


X 


10+03 


1.34 


X 


lO^ 


-03 


1.87 


X 


10+03 


5 


7.21 


X 


10+03 


1.52 


X 


10+04 


1.62 


X 


10+04 


1.85 


X 


lO^ 


-04 


2.77 


X 


10+04 


6 


1.87 


X 


10+05 


3.54 


X 


10+05 


3.51 


X 


10+05 


4.06 


X 


lO^ 


-05 


5.90 


X 


10+05 


7 


6.57 


X 


10+06 


1.25 


X 


10+07 


1.08 


X 


10+07 


1.28 


X 


lO^ 


-07 


1.92 


X 


10+07 


8 


2.74 


X 


10+08 


6.17 


X 


10+08 


4.64 


X 


10+08 


5.31 


X 


lO^ 


-08 


8.87 


X 


10+08 


9 


1.75 


X 


10+10 


3.28 


X 


10+10 


2.52 


X 


10+10 


3.12 


X 


lO^ 


-10 


4.70 


X 


10+10 


10 


1.32 


X 


10+12 


2.40 


X 


10+12 


1.69 


X 


10+12 


2.22 


X 


lO^ 


-12 


3.34 


X 


10+12 



Table 4.3 
X^*^' before and after loop l|4.83|l . 



k 


f>(2fe) 

"o 


before 


^(2fe) 


before 


^{2fc) 


"■0 


after 


j^(2k) 


after 





1.23 


X 


10+02 


2.08 


X 


10+02 


4.84 


X 


10+02 


3.65 


X 


10+02 


4.50 


X 


10+02 


1 


7.81 


X 


10+02 


1.50 


X 


10+03 


3.35 


X 


10+03 


2.79 


X 


10+03 


3.65 


X 


10+03 


2 


3.55 


X 


10+03 


7.74 


X 


10+03 


1.56 


X 


10+04 


1.47 


X 


10+04 


2.00 


X 


10+04 


3 


1.73 


X 


10+04 


3.88 


X 


10+04 


6.63 


X 


10+04 


7.22 


X 


10+04 


1.00 


X 


10+05 


4 


1.85 


X 


10+05 


3.76 


X 


10+05 


4.38 


X 


10+05 


5.37 


X 


10+05 


7.68 


X 


10+05 


5 


3.97 


X 


10+06 


8.59 


X 


10+06 


8.12 


X 


10+06 


9.04 


X 


10+06 


1.40 


X 


10+07 


6 


1.45 


X 


10+08 


2.66 


X 


10+08 


2.46 


X 


10+08 


2.81 


X 


10+08 


4.08 


X 


10+08 


7 


6.62 


X 


10+09 


1.26 


X 


10+10 


1.00 


X 


10+10 


1.19 


X 


10+10 


1.81 


X 


10+10 


8 


3.50 


X 


10+11 


8.00 


X 


10+11 


5.50 


X 


10+11 


6.36 


X 


10+11 


1.09 


X 


10+12 


9 


2.81 


X 


10+13 


5.24 


X 


10+13 


3.81 


X 


10+13 


4.76 


X 


10+13 


7.22 


X 


10+13 


10 


2.59 


X 


10+15 


4.70 


X 


10+15 


3.12 


X 


10+15 


4.18 


X 


10+15 


6.31 


X 


10+15 



We now define 
(4.86) 



pk 



(2fc-2) 




(2k) 



7^(2k-2) 15 -2k 
-"■1 + — ' Pk-1 



k = 0, 



fc > 1, 



so that dHH), glU), and ^AM^ imply 



2k) 
err 



_2fe+2 



2,e 



(2k) 

XX 



(4.87) 



< Vh{\Vo\ + \Vi\) 



Pk' 



-2k-2 



rfe V /i V 16 3 



rk. 



2k+2 



To simplify this expression, we define 



(4.88) 



,,,,, /85 20, 
4 = l5pr^/- + -fc 



and summarize the main result of this section as a theorem. 
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Table 4.4 
Pfc and . 



k 




Pk 









k 










197 


1 


34 


X 


10^ 


-00 


1 





210 


1 


01 


X 


10" 


01 


2 





252 


1 


67 


X 


10" 


02 


3 





288 


3 


58 


X 


10" 


03 


4 





313 


7 


73 


X 


10- 


04 


5 





319 


1 


03 


X 


10- 


04 


6 





305 


5 


96 


X 


10- 


06 


7 





286 


2 


15 


X 


10- 


07 


8 





266 


5 


10 


X 


10- 


09 


9 





248 


9 


15 


X 


10- 


11 


10 





232 


1 


43 


X 


10- 


12 


11 





218 


1 


69 


X 


10- 


14 


12 





204 


1 


58 


X 


10- 


16 


13 





193 


1 


42 


X 


10- 


18 


14 





183 


1 


04 


X 


10- 


20 


15 





173 


5 


98 


X 


10- 


23 


16 





164 


3 


46 


X 


10- 


25 


17 





157 


1 


75 


X 


10- 


27 


18 





149 


6 


86 


X 


10- 


30 


19 





143 


2 


72 


X 


10- 


32 


20 





137 


1 


02 


X 


10- 


34 


21 





131 


2 


94 


X 


10- 


37 


22 





126 


8 


36 


X 


10- 


40 


23 





122 


2 


41 


X 


10- 


42 


24 





117 


5 


40 


X 


10- 


45 


25 





113 


1 


15 


X 


10- 


47 



Theorem 4.8. Suppose k>0, he C'^''+^'^{T), < h{x) <lforO<x<l, and 



£ !i ^o/3. Then the truncation errors V'er? and QerV in (|4.37|) satisfy the hound 



J2k) 



V3 



Q 



(2fe) 
err 



< 



{2k) 
err 



< 



2,E 



{2k) 
err 



-2fe+2 



(4.89) 



< 



2,e 

hm + \vi\) 



{2k) 




PkTk 



2k+2 



, (H3TD . Km . and g: 



where /,„, rk, Pk, and 6k were defined in 

Remark 4.9. The constants in this estimate have been organized to be either 
(1) given in the problem statement or easily computable from h] or (2) difficult to 
compute but universal (independent of h). The first 26 constants in the latter category 
{pk and 9k) are given in Table l44l We have therefore identified the features of h that 
are most likely to affect the validity of the lubrication approximation. In particular, 
higher derivatives are allowed to be large in regions where h is small (since depends 
on the uniform norms of the products j^h^^^d^h rather than on d^h alone). 



Remark 4.10. The term p^Li in the definition of pk ensures that Pj, is a 
nonincreasing function of k. This assumption is useful in the following section when 
deriving a bound on the truncation error of the pressure. Note that pk itself is allowed 
to increase as long as p'^k^'^ does not. It is probably not necessary to include this term 
in the definition of pk since it is not the argmax for 1 < fc < 25, and by that point 9k 
(and hence p'^^'^) appears to be decreasing rapidly without it; see Figure liH 

4.5. Velocity, vorticity, and pressure. We now show how to use the error 
bound we have obtained for the stream function to bound the error in the velocity. 
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p, and log,„6, versus k 



'of 6 









• 

• 


• 

■ ■ • 

• 

• 


• 

• ■ 

• 

• 


• 

• 

• 

• 
• 

• 

• 
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• 

• 

• 

• 


''J - 

• 

• 

• 



10 15 20 
k 



25 



10 


-10 

-20^° 
o 

-30 
-40 
-50 



Fig. 4.1. Plot of and logj^Q 0^ versus k. Note that p^^ initially decreases hut eventually 
grovis almost linearly, indicating that the lubrication expansion is probably an asymptotic series 
rather than a convergent series. The term involving 8k in l|4.89|l is only important when k is small, 
since 8^ converges rapidly to zero as fc — > oo. 



vorticity, and pressure. We define u'^^''\ v^^'^\ uj'^'^^\ and p^^'^) in terms of -0(2'=) as in 
(13. 8|) and define, e.g., 



(4.90) = u:e,act - L^^^''^ 



approxi 



From p.Sp . we tlien liave 

(4.91) (ug^\v^,';^)=Vx^^^4\ 

(4.92) d,p^^J;^=-dyJl^\ 



dyP\ 



(2k) _ _ 
err 




(2fc) _ , 




err 





fc = 0, 

fc > 1, 



wlriclr immediately gives bounds on the error in velocity and vorticity: 
(4.93) 



(2k) 



+ 



£V 



(2k) 



2 \l/2 

.J ^ 


^(2fc) 
'r err 


2,6 




<(•), 


fc > 0, 


^ err 


< 




J,(2fe) 
^ err 


+ £2^+2 

2,e 


TXX 


„ s (•)■ 


fc > 0, 



where (*) is the right-hand side of ()4.89p . Obtaining a bound on the error in the 
pressure is somewhat more difficult, as it relies on the fact that the gradient is an 
isomorphism from L'^(Q,) (the space of square integrable functions with zero mean) 
onto the polar set 



(4.94) 



V° = {f e H^'^inf : (f, u) = whenever u € V} 



where V = {u e i?o(^)^ • V • u = 0}. Given f e V°, there is a unique p G L^{^) 
such that Vp = f ; moreover, p satisfies 



(4.95) 



lbllo<r'|f|-i, /3 = infsup 



l(P,V.u)| 



P u ||p||o|U|i 
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Here we use a standard (unweighted) Sobolev norm for f. More precisely, as || ■ ||i 
and I'll are equivalent on HQ{fl), the negative norms 

(4.96) |f|_i= sup |(f,(u;«))|, ||f||_i = sup |(f,(u;«))| 

|«|? + |«|2 = 1 ||„||2 + ||„||2^1 

are equivalent on H~^{n)'^. Since |u|i < for all u E ifo(r2), we have |lf|l-i < 

|f|_i for all f € H-\n)^. 

Explicit estimates [71 [23l [10] for the LBB constant {3 in (|4.95p have been obtained 
for rectangular domains (with no periodicity), e.g., 

(4.97) isinJ</3<^, £ = max(^^,^), 17 = (0, Li) x (0, L2). 

The lower bound here also works for an x-periodic rectangle as the condition that 
u|a;=o = and u|a;=ij = is more restrictive than u|a;^o = '^\x=Li- Explicit estimates 
are also known for domains that are star-shaped with respect to each point in a ball 
of radius R contained inside fi; see [T31 [13] • Our interest in the present work is in 
x-periodic domains with the upper boundary given by a function h{x). Such domains 
are not in general star-shaped, so the previously known results do not apply. In 
[24j . we improve the estimate for the lower bound on /3 in (|4.97|) for an a:-periodic 
rectangle by a factor of about 3.5 and show how to avoid invoking Rellich's theorem 
in the change of variables to the case that il is x-periodic with the upper boundary 
given by h{x). It is shown that 

<-(l+r-^) fh.\ "max (^4——^ /lo - mino<,<L /i(a;), 

^ -i^^^'^>[ho) ^['ho'hoj' /ii=maxo<.<iM^), 

where tq = max(||ft,j;||oo, |i ^/i/iajxlloo^)^"'^ and L is the period of h(x). In the current 
case, the length scales H and W were chosen so that L = 1 and hi < 1 in the 
dimensionless problem. Thus, solving Vp = f yields 



(4.98) ||p||o<rlf|-i, <max(^9V ' j (l + '^i 

The dependence on gap thickness /iq occurs because p can change rapidly in the gap 
without a large penalty from u. In [24], an example is given to show that the factor 
of h^^^'^ in the formula (|4.98l) for cannot be improved. 

(2k) 

We have reduced the problem of bounding perr to that of bounding the functional 
on the right-hand side of Vpir? — ffc in (|4.92p . namely, 

(499) j,^.^^^ _ f/nwl?Uy ~e^a;ea;actWa;dA, fc = 0, 

First, we check that belongs to . If w, u G HQ{Vt) and Ux + Vy = 0, the function 
(/)(a;, y) = u{x, rf) drj satisfies (j)y — u, (f>x — —v and so belongs to ^. As a result. 



(fo,(w;u)) = I -UJ^^^cjyyy +UJ exact {(f'vv+e'^^'xx) dA 

(4.100) 

^yj'l^yy + (A.i: exact) iAe<j)) dA = 0, 
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and for fc > 1, 



Next, we bound the norm of f^. From (|4.99p . we see that 



|fol-i < 



,(0) 



^ ^ exact n 1 



< 



(2k} 



+ 



omit if A; — 1 



£2 (2fe-2) 



fc > 1. 



Denoting the right-hand side of (|4.89p by (*), we claim that 



(4.101) 



6 £ 
^|lWe2;acf|lo < (*), (^ = 0), 
Tn rt 



<W, (fc>i)- 



'0 'fc 

Once this is shown to be true, we will have the following bound on |ffe|-i: 

(4.102) |ffe|_i < (l + r2)(*), (fc>0). 

Together with (14.981) and the fact that < tq for fc > 0, this will give 



(4.103) 



p(2fc) 



< max 



(9v^^^^o^^^) (^.+^fc-^)'(*)- 



Note that there is at least a power of r^T^ in (*) to prevent this bound from diverging 
as Tfe — ^ cx). It may be possible to improve the bound in this regime by replacing 
||e^Wej;(jct|lo in (|4.101|) by We'^dxUJexactW-i, but this seems very difficult. At any rate, 
if Tfe = oo, then (*) = 0, h(x) is a constant function, the exact and approximate 

vorticity are constants, fk is the zero functional, and Pexact — P^err = 0. Let us now 
prove (|4.10ip . For fc > 1, this follows from 



(4.104) 



Pk- 



-2k 



15 



h 85 20, 



rk-1 V A V 16 3 



rk-i 



2k 



< 



Pk 



-2k-2 



Tk V /i V 16 3 



2k+2 



which holds because p^^"^ and are nonincreasing functions of fc; see Remark 14.101 
and the definition of in (|4.6ip . For fc = 0, we use Theorem 14.41 to conclude that 



(4.105) 



90 



1/2, e 



91 



1/2, £ 



where 50(2^) = Vo and 51(2;) = (1 + e'^hx{x)'^) Now 



90 



(4.106) 



1/2, e 
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and 
(4.107) 



91 



1/2, e 



2u \-l/2 



= vn {h-') {■)-' + 



dx 



2 / 2 ^^xx 



dx 



(4.108) 



1 4 \ / \ 2 



^1+1 2+^4 



Since we have assumed that e < 7'o/3, we conclude that 



(4.109) ^WiOexactWo < V h{\Vo\ + \Vi 



Comparmg this to (|4.87p with k = and noting from Table IT4l that > 15, we 
obtain ()4.10ip as claimed. Thus, we have proved the fohowing theorem. 

Theorem 4.11. Suppose k>0,h€ C'^''+^''^(T), < ho < h{x) <l for x eT, 
and e < ro/3. Then the truncation errors of the stream function, flux, velocity, 
vorticity, and pressure satisfy the bounds 



(2k) 

err 



(2k) 



ev 



(2k) 



1/2 



(4.110) 



I2k) 



„(2fe) 
F err 



< max 



(9V^/^^V^/^) irk+r-^Yi^), 



where T — [0, l]p is the periodic unit interval 
(4.111) 

(4.112) rfc 



= V^id^ol + I^il) 
1 



max 

l<f<2A;+2 



h'~'dih 



e_ h 

' Tk V h 



e 



PkTk 
dm 



2k+2 



h{x)-"^ dx, 



and pk. Ok are constants independent of h that can he computed once and for all as 
described in section W^ and listed in Table 

5. Finite element validation. In this section, to test the error bounds of The- 
orem |4TT1 we compute Werr\ Verr^-) P^err\ UJerr^ numerically for the simple geometry 
described by 



(5.1) 



h{x) 



1 + a 1 — a 



sin(27rx). 



Case 1: a = 1/5, 
Case 2: a =1/100, 



(2k) 



with boundary conditions Vq — —0.5, Vi — 1. We do this by comparing Uapprox, 
v'iijl^Jrox, P^approx, ijj'approx in (|4.90p to finite element solutions of the Stokes equations 
on appropriately rescaled geometries. 
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Actual errors 




8 



Error bounds from Theorem 4.1 1 




£ 

Fig. 5.1. Top: plot of ||ul.r?' = (||«l?r ' |jf + Wev'-^r ^ Wl^^y^^ , ||pe?r'||o, and ||w^r?''||o f"i~ 
a = 1/5, 0.04 < e < 0.3, and 2k £ {0, 2, 4, 10, 20}. The slopes of the lines were computed via linear 
regression using the smallest 10 values of e for which the finite element solution is trusted (e > 0.066 
for 2k = 10 and e > 0.09 for 2k = 20). As expected, for fixed k, the error is 0(£2fe+2). Bottom: plot 
of the error bound (*) in l|4.11ip . using Vo = -.5, Vi = I, Ii = 2.236, I3 = 24.60, and rt = 0.3559 
for k >0, as appropriate for h{x) in i5.lt with a = 1/5. 

The results are sunimarized in Figures I5.1H5.5I For 21 values of e spaced expo- 
nentially between Eq = 0.04 and £20 = 0.3, we set up a logically rectangular, M x N 
finite element mesh on the domain 



(5.2) 



n,^{{x,y) : 0<a;<l, < y < eh{x)}. 
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po 

0.15 

0.1375 

0.125 

0.1125 

0.1 

0.0875 

0.075 

0.0625 

0.05 

0.0375 

0.025 

0.0125 







p4 






0.00018 






0.00015 






0.00012 






9E-05 






6E-05 






3E-05 













-3E-05 






-6E-05 






-9E-05 






-0.00012 




1 


-0.00015 






-0.00018 





plO 






1 .2E-06 






1E-06 






8E-07 






6E-07 






4E-07 




2E-07 






1 


-2E-07 




-4E-07 




-6E-07 



p20 
1 .5E-07 
1E-07 
5E-08 


-5E-08 
-1E-07 
-1.5E-07 
-2E-07 



Fig. 5.3. Contour plots of pexact, Pern Perr, Perr , o,nd Perr foT h(x) in i5. ID with a = 1/5 and 
e = 0.099. The "exact" solution was computed using a least squares finite element method with 15 
node quartic triangular elements on a 2208 x 96 grid. 
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0.04 0.05 0.07 0.1 



0.2 0.3 



Fig. 5.4. Error estimates and actual errors with a = 1/100. 





1 

0. 
0.6 
0.4 
0.2 




colO 

— I 1.20E-05 
9.00E-06 
6.00E-06 
3.00E-06 

IO.OOE+00 
-3.00E-06 
-6.00E-06 
-9.00E-06 
-1.20E-05 



plO 

' ' -06 
■06 
■06 
■07 





2.50E- 




1 .94E- 




1.39E- 




8.33E- 




2.78E- 




-2.78E- 




-8.33E- 


i 


-1.39E- 




-1.94E- 




-2.50E- 



Fig. 5.5. Plots of uj exact, Pexact, and p^JrV with a = 1/100, N - 



M = 4800. 
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The mesh points are ahgned verticahy with equal spacing Ay — h{x)/N, while the 
grid spacing in the x-direction is chosen to keep the aspect ratios of the grid cells as 
close to 1 as possible; we do this by solving an ODE to enforce Aa; sa h{x)/N, which 
also determines M. For a = 1/5, we use iV = 96 with M ranging from 768 to 5376 
as e ranges from 0.3 to 0.04; for a — 1/100, we use iV = 64 with M ranging from 
1600 to 10368. Four-by-four blocks of neighboring grid cells are merged and cut into 
two 15 node triangles. Interior nodes of the triangles are adjusted to keep the edges 
straight except on the top boundary, where we use quartic isoparametric elements. 
We solve the Stokes equations on this mesh using a least squares finite element method 
similar to j6] but using quartic elements to model the velocity components u and v, 
the pressure p, the vorticity uj = — Uy, and two strain rates t = Uy + and 
^ = Vy — Ux- We use multigrid to solve the resulting system of equations, which takes 
from 3 to 15 minutes on a 2.4 GHz desktop machine with 16 GB RAM. 

Once the finite element solution is known at the grid points, we normalize the 
velocity, pressure, and vorticity as described in section [2] and rescale the domain 
from to f2. We then use the method described in Appendix |X] to compute 
V'^^-', ■ . ■ , V'*-^"'' and their derivatives through order 3 at the grid points. Next, we use 

and p^"^^^ for = 0, . . . , 10. For 



the formulas in (13.81) to obtain u^'^^\ v^'^^\ uj^'^'^ 



(2k) 

pressure, we use 20 point Gaussian quadrature to integrate px along the x-axis to 



determine p'^'"'^ {x,0) at the mesh nodes. The integration oip^y'^^ in the y-direction is 
done analytically. With the expansion coefficients in hand, we evaluate 



(5.3) 



,(2fc) _ 



exact ^ 



(2k) 
approx") 



approx 



,,(2fc) 

err 


2 


£„(2/c) 


2 

l,e 


„(2fc) 
r err 


2 

7 




^(2/c) 

err 



etc., at the grid nodes, where we use the finite element solution for u^xact- We then 
run through the triangles and sum up the local contributions to the errors 



(5.4) 



by interpolating the values at the grid nodes and integrating the resulting polynomials 
on the triangle; this step is very similar to the assembly of the stiffness matrix. Finally, 
we store the results in a file for visualization (see Figures 15.21 15.31 and 15. 5p and 
record the norms of the truncation errors for comparison with the error bounds of 
Theorem gm 

The results of this comparison are shown in Figures 15.11 and 15.41 As expected, 
for fixed fc, the actual errors decay as 0{£^^^'^). The a priori error bounds eventually 
decrease like 0(e^'^"'"^) as well, but the term involving 6^ in (|4.11ip is significant over 
this range of e in some of the cases, causing the slopes to be larger: 





fc = 


fc = 1 


fc = 2 


fc ^ 5 


fc = 10 


a = 1/5 


12.5 


0.94 


0.16 


.00096 


1.3 X 10-" 


a = 1/100 


257 


19.3 


3.2 


0.020 


2.7 X 10^1° 



Ok_ 



This effect is much more pronounced when a = 1/100 in (|5.1I) due to 

fh_^l /313^ [3.32, a = 1/5, 
^ ' V/i 2V2^a^2a2 |61.4, a = 1/100. 

The deviation from linearity in the plots of "actual error" for small e and large fc is 
due to error in the finite element solutions, which are accurate to about 9 digits. This 
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occurs sooner when a — 1/100 since the pressure and vorticity of the exact solution in 
the vicinity of the narrow gap increases as a decreases, and also because we were forced 
to use a coarser mesh with a = 1/100 to avoid running out of computer memory in the 
finite element simulations. The data points with e = 0.099 in Figure 15.11 correspond 
to the contour plots in Figures 15.21 and 15.31 where we plot Uexacu ^^err\ Pexacti and 
Perr'' for 2fc = 0, 4, 10, 20. The data points with e — 0.099 in Figure [Ol correspond to 
the contour plots in Figure [531 We remark that the apparently large value of pir? in 
the narrow gap in Figure 15.51 is due to smoothing in the least squares finite element 
solver; the expansion solution is more accurate than the finite element solution in this 
region of the domain. The error patterns that emerge in all these cases are rather 
interesting, indicating that the spaces 'H2k in Theorem 13.31 (the structure theorem) 
can be quite complicated even for simple curves h(x). 

Although our estimates for the error in pressure include an additional factor of 
/iq '^''^(/'fe all our numerical experiments (including complicated geometries in 

3/2 f 2 ^ 

which the inf-sup constant /?^^ does exhibit /ig behavior) indicate that ||perr |lo is 
comparable to Hwir? ||o- In fact, for large fc, pressure seems to be the most accurately 
computed variable; see Figures 15.11 and 15.41 We do not know how to explain this as 
the pressure is determined by solving (|4.92p . which involves inverting the operator 
V : L'^iyi) -> i/~^(Sl)^. For some reason, in lubrication-type problems, the right- 
hand side ffc belongs to a subspace of H~'^{VlY that is not amplified by (V)~^ when 

(2k') 

solving Vperr = f fe . 

The following table shows the minimum ratio of the a priori error estimate to the 
actual error Ijuir? || i.£ for the data points in Figures [5.11 and [5.41 that were used to 
compute the slopes of the best-fit lines: 



k = 





1 


2 


5 


10 


(min ratio, a = i/5)^/(^'^"'"^) 


11.1 


7.0 


5.0 


2.8 


2.3 


(min ratio, a = 1/100)^^^^''+^' 


34.0 


8.6 


5.5 


3.0 





For example, in the 10 calculations (with e ranging from 0.04 < e < 0.099) that were 
used to determine the slope of the 2k = A line in Figure 15.11 the ratios of the a priori 
errors to the exact errors rang ed between 1.608 x 10^ and 1.617 x 10'', so we recorded 
v^l.608 X 10^ « 5.0. This table gives information on how far the values pk in Table 
are from their optimal values. For example, if we increased by more than a factor 
of 2.8 while holding 65 fixed, the estimate (j4.110p would fail to hold for this geometry. 
Since r^^ in ()4.61|) is used as a convenient upper bound on all the integrals l-E^^^jl^/^^ 
and \eI^^^\^/^^ that arise in the definition of Q'^^'^^ and also in the bounds for ||'0iaf''||o 

and Ij/i^V'i^^EK^'' llo: it is remarkable that the values of pk we computed are within a 
factor of 3 of optimal for A: = 5, A: = 10, and perhaps all fc > 5. 

6. Discussion. Although we are able to estimate the effective radius of conver- 
gence pfcTfc quite closely, our estimates of HV'err'' lU.e, H'^irr'' |jo, etc., are likely to be 
several orders of magnitude too large. One shouldn't expect an a priori bound that 
holds for all geometries alike to provide an exceptionally sharp bound for any spe- 
cific geometry. Instead, our analysis provides a clear picture of the features of h{x) 
that cause the effective radii of curvature rkpk to become small, namely, large values 
of h^~^d^h. No previous study has ever described how the constant hidden in the 
Q^^2k+2-^ depends on h; instead, h has always been fixed at the outset and only the 
limit as e — > has been considered. 
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Approximation of velocity at a point 




(2k) ^ ^ 

Fig. 6.1. Comparison of u approx (solid lines) to Uexact (dots) at the point {x,y) = (g, j/i) for 

2k = 0, 2, 4, 6, 10, 16, 20, 30, 50. Here h{x) = | + | sin(27ra;), Vo = -0.5, and Vi = 1. This function 

h(x) is real analytic and periodic, yet the expansion solution appears to be an asymptotic series 

rather than a convergent series. 

Another feature of this analysis is that it separates the constants into two types: 
those that are (1) given in the problem statement or easily computable from h\ or 
(2) difficult to compute but universal (independent of h). We listed the first several 
constants in the latter category {pk and Ok) in Table l44l It is interesting that pk 
actually increases until 2k = 10 and doesn't get as bad as po again until 2k = 26. 
However, at that point it seems to be decreasing steadily like 1/k, indicating that 
the effective radius of curvature in our a priori error bound will shrink to zero as 
A: — >■ oo. The reason for this is that the recurrences (|A.3|) and (|A.4p relating the 
matrices and S^^*^) to their lower order counterparts cause the norms of these 

matrices to grow like kl. Thus, although pk involves kth roots of these constants, 
these fcth roots still grow linearly in k. On the other hand, if h is real analytic as well 
as periodic, a standard contour integral argument shows that there is an r > such 
that ||9^/i||oo < k\ for all A; > 0; thus, the constants will remain bounded away 
from zero. For example, if h{x) is of the form (jS.ip . one may show that if a S (0, 2/3], 
then the largest value of \\j\h^~^d^h\\\lf' occurs when ^ = 2, so all the are equal 
to tq = (tt-^/I — a)^^ . It is conceivable that when h is real analytic, the norms of the 
functions ■0^^'^^ grow slowly enough that the stream function expansion converges in 
spite of the fact that the matrices Af''^'^ and B^'^'^^ in their representation p.22p blow 
up like fc!. This would simply mean that we chose a bad basis in terms of which to 
represent ?/;. We used orthogonal polynomials in (j4.77l) to improve this basis, but there 
may be other improvements. Figure |6T] shows that this is not the case. Even when h 
varies sinusoidally, the expansion solution appears to be an asymptotic series rather 
than a convergent series: all the variables, including the flux terms Q^'^^\ appear to 
grow like fc! as k becomes large. 

Nevertheless, the expansion solutions can be extremely accurate (almost exact) as 
long as they are used for a geometry that falls within the effective radius of convergence 
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of the truncated series. It is hoped that the estimates in this paper wiU help to identify 
these cases and provide practical a priori (as well as a posteriori) error estimates for 
many interesting problems. 

Appendix A. Implementation. We have developed two methods for comput- 
ing the higher order corrections described in section using a computer. In the first, 
we use Mathematica to evaluate the derivatives and antiderivatives in recursion p.6p 
and Algorithm 13.11 symbolically. With this approach, the main challenge occurs at 
the step where Q^^*^") is defined as a definite integral. We do this through pattern 
matching and symbol replacement. At the stage where the definite integral is to be 
evaluated, we replace all instances of d^h in the integrand by j! tj/h^~^. Each term in 
the result (call it R) will contain a factor of or ft,^^, with no other dependence on 

(2k) 

h. For each k — ko, . . . ,0 and j = 1, . . . ,d2k, we find the terms in R that contain (pj 
(left in the form ■ ■ ■ described in Algorithm 13.21) as a factor. These terms are 
removed from R while their symbolic integrals (with 93^^'^^ //i™ replaced by ImE^'^j) 
are divided by 2J3 and added to the desired flux Q^'^^°\ By running through the 
in decreasing k order, we convert higher order products (e.g., t\t2/h?) into symbols 
(e.g., I^e'^^I) before one of their lower order factors can be converted incorrectly (e.g., 

into /a-Es 1^2)- This approach is effective through 6th or 8th order but becomes rather 
slow as the complexity of the expansion increases. 

The second approach is much faster and can be implemented in any modern pro- 
gramming language. We have written a version in C'^^ and a version in Mathematica. 

(k) 

Instead of representing the basis functions t/?^ for Tlk using a computer algebra sys- 
tem, we represent them as [k + l)-tuples of integers. For example, the functions 1, 
hx, \h?hl.hxxx, and j^h'^hxhxxhxxxx in "Ho, 'Hi, "He, and 'Hr are represented by (0), 
(0, 1), (2, 3, 0, 1, 0, 0, 0), and (4, 1, 1, 0, 1, 0, 0, 0). A tuple (iq, . . . ,ik) represents a basis 
function for Hfe iff 

(A.l) ii+2i2-\ \-kik^k, = ^2 + 2^3 H h (fc - l)u.- 

We begin by constructing the basis sets $fc for < fc < 2fco and storing them as 
(fc -|- 1) X dk integer matrices with columns corresponding to the ^j''^- This is done 
using Algorithm [321 which returns the columns sorted lexicographically from the last 
slot to the first slot (e.g., (3,0,3,0)"^ < (2,3,0,1)^ < (3,1,1,1)'^). Sorted columns 
allow us to find the column index corresponding to a given tuple in log2 dk time. 

Next, for < fc < 2fco — 1, we compute the operators hdx and hx- from Hk to 
Hk+i and store them as sparse integer matrices of dimension dk+i ^ dk- If column J 
of $fc contains the tuple (io, . . . , ik), we define ik+i = and compute 

(A. 2) hx- : (io, . . . , ifc) H> (io, ii -I- 1, 12, . . . , ifc+i), 

hdx : {io,- . ■ ,ik) '-^ ^ v(r- + l)(io -I- 1, . . . , v - 1, V+i + 1, ■ • • , ifc+i), 

where the omitted indices are unmodified and the -1-1 and —1 cancel in the first slot 
when r = in the sum. The factor of (r + l) is due to the factorials in the definition of 
the </5j'^''- The column index I of each (fc + 2)-tuple in the result is found in <I>fc+i, and 
the corresponding coefficient (1 or ^(r + 1)) is added to the Ith. row and Jth column 
of the sparse matrix representing hdx or hx-. The entries of these sparse matrices 
are positive, and the column sums (i.e., 1-norms) are all equal to 1 for hx- and to 
io -I- 2ii H h (fc + l)ik = 2fc for hdx (by (|A.1|) V 
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Once the operators hd^ and h^- are known, we use them to recursively compute 
the matrices A^^^) = VqA^^'''^ + ViA[^''^ and B^^^) in (jX^ . We start by setting 
A^°^ = (0, 1, -2, 1)"^, Af^ = (0, 0, -1, 1)^, and St") = (0, 0, 3, -2)^ as in Example^ 
For 1 < fc < fco, we mimic the proof of Theorem 13.31 to build up A'^^''^ and S'^fc) ^ow 
by row. For 4 < n < 2k + 3 and i = 0, 1, we use sparse matrix-vector multiplication 
to define the rows 

, -2[hd,~{n~2)hShd,~in~3)h,]\A^''-^\n,:) 

A^'Hn,:) = 

(A.3) 

^(2fe)^^ ,^^f -2[hd.. -in- l)h.,][hd.,, -in- 2)K] :) 

' \ "(^ ~ 1) / 

If fc > 2, then for 6 < < 2A: + 3, we add the following vectors to A[^''\n, :) and 
B^'^^\n,:), respectively: 



n(n — 1) 



, T 
Tl \ 



(A.4) 



-[ha^-(n-2)ha,][hd^-{n-?,)h^][hd,:-(n-i)h^][hd^-{n-b)h^] [/if''"*' (n,:)^] ' ^ 
n(n-l)(n-2)(n-3) 

-[/ia^-(«-l)/i^][/i9^-(n-2)/i^][feO^-(n-3)/i^l[/t3,;-(K-4)/i^][B<^''^^>(».:)^ ' 
n(n-l)(n-2)(n-3) 

(2fc) .(2fc) 



Next we zero out rows and 1 of A^ , A\^ , B^'^^\ and set 

2fc+3 2fc+3 

Ar)(2,:)= 5:(n~3)Ar)(n,:), (3, :) = J] (2 - :), 

(A.5) 

^(2fe) :) ^ ^ „ 3)ij(2fc) ^n, :), i?*^^) (3, :) ^ ^ (2 ~ n)B^^''^ in, :). 

n=4 n=4 

Finally, we subtract ( J,^^) from ^^^'^•'(2, 1) and add it to ^^^'^-'(3, 1) to account for 
the boundary data, where we recall that the rows and columns are indexed starting at 
and 1, respectively. Using this approach, our C++ code can compute these matrices 
through order 2k — 50 using floating point arithmetic in a few seconds, while our 
Mathematica code can compute through order 2k — 30 in exact rational arithmetic 
in about an hour. This allows us to explore the properties of the stream function 
expansion and test our error estimates to quite a high order. 
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